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Abstract 

This  study  identified  techniques  and  software  available  for  the  optimization  of  doubly  curved 
shells  and  applied  them  in  the  context  of  a  large  nozzle  shape.  An  optimality  criteria  scheme 
that  can  reduce  solution  time  was  evaluated  and  compared  to  the  Method  of  Feasible  Directions. 
MSC/NASTRAN  and  ASTROS  were  used  to  perform  finite  element  analysis  and  optimization,  and 
the  results  were  compared  to  theory.  The  programs  give  virtually  identical  results,  and  if  plates 
and  shells  are  carefully  modeled,  then  stresses,  displacements  and  modes  are  accurate  to  within  ten 
percent.  A  Mindlin-type  axisymmetric  finite  element  was  implemented  in  ASTROS  that  preserved 
accuracy  and  reduced  the  size  of  the  stiffness  matrix  by  a  factor  of  four. 

Nozzle  optimization  was  performed  using  static  pressure  and  thermal  loads,  constrained  by 
the  Von  Mises  stress  criteria.  Software  errors  in  ASTROS  were  documented,  and  four  characteristic 
stress  regions  identified  for  the  optimized  nozzle. 
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A  COMPARISON  OF  THE  OPTIMIZATION  AND 


ANALYSIS  OF  DOUBLY  CURVED  SHELLS 
USING  MSC/NASTRAN  AND  ASTROS 


I.  INTRODUCTION 


Purpose  And  Objectives 

Rigorous  optimization  of  aerospace  systems  is  required  to  extend  range,  reduce  costs,  and  meet 
preformance  objectives.  The  trend  toward  complex,  high  performance  systems  requires  early  sys¬ 
tem  level  multidisciplinary  optimization.  This  has  driven  the  implementation  of  multidisciplinary 
optimization  schemes  in  finite  element  analysis  software. 

Doubly  curved  shell  structures  are  a  key  element  in  aerospace  systems.  Their  high  strength 
has  led  to  use  in  aircraft  skins,  rocket  fuel  tanks  and  nozzles,  among  other  things.  The  purpose  of 
this  study  is  to  identify  techniques  available  for  the  optimization  of  doubly  curved  shells  and  apply 
them  to  an  aerospace  structure.  This  will  be  done  by  researching  robust  optimization  algorithms 
and  the  software  that  implements  them.  The  software  will  be  tested  for  accuracy,  then  used  to 
optimi  •  the  wall  thickness  of  a  large  nozzle. 

Literature  Search 

Literature  In  Optimization:  Structural  weight  optimization  is  important  to  the  design  of 
Aerospace  vehicles  since  weight  has  a  great  impact  on  vehicle  performance.  Optimization  is  gen¬ 
erally  performed  in  the  context  of  constraints,  and  various  methods  have  been  developed  to  deal 
with  them.  The  constrained  optimization  problem  may  be  changed  to  an  unconstrained  problem 
by  transforming  the  cost  and  constraint  functions  into  a  new  cost  function  that  penalizes  constraint 
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violations.  The  problem  is  then  optimized  using  unconstrained  techniques.  Arora  (2)  details  tech¬ 
niques  that  fall  into  this  category,  and  their  inherent  instabilities  near  constraint  boundaries.  More 
popular  and  powerful  methods  use  constraints  directly  in  the  optimization  process  using  Lagrangian 
multipliers.  Miura  (17)  describes  a  well  conditioned  Modified  Feasible  Direction  algorithm  that  is 
faster  than  steepest  descent  methods.  Venkayya  (26,  27)  details  an  Optimality  Criteria  Method  that 
resizes  or  scales  the  design  based  on  Kuhn-Tucker  conditions  and  the  location  of  the  constraints. 

Literature  In  Plate  And  Shell  Structures:  In  general,  a  plate  structure  is  defined  as  the  solid 
material  enclosed  between  two  closely  spaced  planer  surfaces,  and  a  shell  structure  as  the  material 
between  two  curved  surfaces  (7:Gibson).  The  utility  of  plate  and  shell  structures  is  illustrated  in 
creation.  Insect  exoskeletons  are  commonly  made  of  thin  shells,  and  plate  like  structures  are  seen 
in  the  plant  world.  Plate  and  shell  structures  have  a  high  strength  to  weight  ratio  that  is  useful  in 
aerospace  applications.  Nozzles  are  one  important  application  of  doubly  curved  shells. 

Closed  form  solutions  to  the  structure  and  vibration  problem  exist  for  a  number  of  simple 
plate  and  shell  structures.  Timoshenko  (23)  has  solved  for  stresses  and  displacements  in  a  uni¬ 
formly  loaded  round  plate  with  uniform  thickness  and  clamped  edges.  Meirovitch  (16)  illustrates 
the  derivation  of  natural  frequencies  and  mode  shapes  for  this  problem.  A  spherical  dome  with 
rigidly  fixed  edges  is  another  doubly  curved  shell  structure  for  which  a  closed  form  solution  exists. 
Timoshenko  (23)  illustrates  both  the  derivation  of  the  exact  solution  and  a  simpler  approximate 
solution  which  was  first  developed  by  Hetenyi.  Gibson  (7)  illustrates  the  FORTRAN  coding  of  the 
approximate  solution. 

Literature  In  Nozzles:  Supersonic  convergent/divergent  nozzles  are  used  in  a  variety  of  ap¬ 
plications  including  rocket  and  turbine  engines  and  supersonic  wind  tunnels.  Nozzle  shape  and 
length  may  be  optimized  based  on  the  characteristics  of  the  flow  and  mission  profile.  Oates  (19) 
illustrates  the  relationship  between  expansion  ratio  and  thrust  at  different  altitudes.  The  nozzle’s 
structural  design  is  driven  by  material  strength  and  design  stability,  subject  to  pressure  loading 
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and  other  mission  loads.  By  assuming  one  dimensional  flow,  the  pressure  loading  may  be  approx¬ 
imated  using  the  Laval  nozzle  equations.  Kuethe  and  Chow  (12)  detail  the  development  of  these 
equations.  Mission  loads  may  include  transient  thermal  effects  from  the  gas  flow,  thermal  effects 
from  nuclear  bursts,  erosive  effects,  and  a  vibration  spectral  density  environment.  For  instance, 
erosion  is  the  primary  concern  in  analyzing  solid  propellant  rocket  nozzles,  followed  by  thermal 
stresses  and  conductivity  (6:Galati).  Th"  steady  state  operating  temperature  of  the  nozzle  is  an 
important  consideration  for  liquid  propellant  rocket  engines.  Barrere,  Jaumotte,  DeVeubeke  and 
Vandenkerckhove  (3)  illustrate  methods  developed  by  Ciniaref  and  Dobrovolski  for  calculating  these 
temperatures  in  a  liquid  propellant  rocket  nozzle  with  regenerative  cooling. 

Literature  In  Finite  Elements  And  Software:  With  increasing  demands  for  vehicle  perfor¬ 
mance  there  arose  the  need  for  multidisciplinary  optimization  early  in  the  design  process.  In  the 
past,  large  software  packages  performed  analysis  in  one  discipline  at  a  time.  It  was  up  to  the 
user  to  collect  the  necessary  design  sensitivity  information  from  the  various  disciplines  and  put 
them  in  a  form  that  could  be  optimized.  The  Air  Force  Wright  Aeronautical  Laboratories  initiated 
development  of  ASTROS,  short  for  Automated  Structural  Optimization  System.  This  is  a  large 
software  package  that  collects  static,  modal,  aerodynamic  and  dynamic  analysis  together  for  mul¬ 
tidisciplinary  optimization  (10:Johnson).  The  MacNeal  Schwendler  Corporation  developed  a  new 
MSC/NASTRAN  solution  sequence,  SOL  200,  that  collects  static,  dynamic,  and  buckling  analysis 
together  for  the  same  purpose  (17:Miura). 

Finite  elements  for  shells  have  been  difficult  to  develop.  Curved  elements  closely  model  the 
physical  situation  but  are  complicated  Modeling  the  shell  as  a  set  of  finely  faceted  flat  elements 
is  less  complex  and  works  well  enough  to  compete  with  curved  elements  (4:Cook).  An  alternative 
for  problems  with  axial  symmetry  is  an  axisymmetric  shell  element.  Cook  details  the  development 
and  accuracy  of  a  Mindlin  type  axisymmetric  shell  element  with  membrane,  transverse  shear  and 
bending  stiffness.  MSC/NASTRAN  (14:MacNeal)  incorporates  a  conical  axisymmetric  shell  ele- 
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ment  with  all  but  drilling  degrees  of  freedom  at  each  nodal  ring.  A  Fourier  analysis  scheme  is  used 
to  obtain  the  response  to  non-axisymmetric  loads. 
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II.  OPTIMIZATION  BACKGROUND 


The  General  Optimization  Problem 

Optimization  has  been  pursued  in  a  variety  of  disciplines,  but  regardless  of  the  discipline 
the  optimization  problem  can  be  reduced  to  a  general  form.  Numerous  mathematical  optimization 
techniques  have  been  developed  based  on  this  standard  form,  which  is  defined  by  Arora  (2)  as 
follows: 

Given  the  design  variables 

*  =  (ii,...,sr„)  (1) 

and  the  function 

F(x)  =  F(xu..  ,x„)  (2) 

find  the  x  that  will  minimize  F(x)  subject  to  the  following  constraints: 


Zj(xi..  ..,xn)<0j=l...k  (3) 

Zj(x  =  0  j  =  k+\...l  (4) 

Xj  Xj  upper  <  0  J  —  1  ...  71  (5) 

Xj  Xj  lower  ^  0  j  —  \  ...  71  (6) 


Equation  3  represents  the  k  inequality  constraints,  Equation  4  the  l  —  k  equality  constraints,  and 
Equations  5  and  6  represent  the  side  constraints  on  the  design  variables.  The  function  F(x)  may 
be  referred  to  as  the  cost  or  objective  function  (2:Arora).  In  general,  the  objective  function  and 
constraints  may  be  nonlinear,  implicit  functions  of  the  design  variables. 

The  design  variables  define  an  n  dimensional  design  space.  Any  arbitrary  xv  is  a  design.  If 
the  x'  produces  inequality  constraints  that  are  more  positive  than  some  user  defined  constant  e,  or 
equality  constraints  that  are  not  zero  to  within  ±e,  the  constraints  are  violated.  Similarly,  inequality 
constraints  between  zero  and  —  e,  and  equality  constraints  between  ±e  are  considered  satisfied  and 


active.  Inequality  constraints  that  are  more  negative  than  —  e  at  the  design  are  characterized 
as  satisfied  and  inactive.  The  activity  or  inactivity  of  a  constraint  relates  to  its  involvement  in 
moving  about  the  design  space.  At  a  given  design,  active  constraints  provide  important  information 
about  the  direction  toward  the  optimum.  Inactive  constraints  remain  part  of  the  problem  but  the 
information  they  provide  is  less  critical.  For  numerical  solutions  the  user  defined  variable  e  is 
necessary  because  in  a  practical  sense  it  is  impossible  to  exactly  satisfy  a  constraint.  Without  some 
finite  e,  all  the  constraints  would  be  characterized  as  violated  or  inactive,  and  the  optimization 
routines  would  have  to  act  on  all  the  variables.  This  would  be  inefficient. 

Given  the  design  a;1',  we  may  characterize  it  as  feasible  or  infeasible  based  on  the  condition 
of  the  constraints.  If  there  are  any  violated  inequality  or  equality  constraints,  and  the  design  does 
not  lie  within  the  bounds  specified  by  the  side  constraints,  then  the  design  is  infeasible.  On  the 
other  hand,  if  all  the  constraints  are  satisfied  then  the  design  is  feasible.  A  feasible  design  may 
have  some  active  constraints,  in  which  case  it  is  a  constrained  optimization  problem.  Otherwise  it 
is  an  unconstrained  problem.  At  any  given  xv  there  are  three  possible  design  states; 

1.  The  design  is  feasible  and  there  are  no  active  constraints. 

2.  The  design  is  feasible  with  one  or  more  active  constraints. 

3.  Thf  design  is  infeasible  since  one  or  more  of  the  constraints  are  violated. 


The  design  state  is  important  in  determining  the  action  to  be  taken  during  the  optimization  process. 

The  intent  of  the  design  optimization  process  is  to  find  the  feasible  design  x*  that  minimizes 
the  objective  function.  To  do  this,  the  conditions  that  must  exist  at  the  minima  need  to  be 
identified.  If  the  minima  is  bounded  and  the  problem  unconstrained,  then  the  familiar  calculus 
concept  of  the  vanishing  gradient  at  the  minima  is  a  necessary  condition,  or: 


VF(x) 


dF(x)  dF(x ) 
dxx  ’  "  dxn 


(7) 
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Arfken  (1)  has  demonstrated  that  if  the  design  at  the  minima  includes  active  constraints,  then 
the  are  no  longer  arbitrary  and  in  general  do  not  vanish.  This  leads  to  formulation  of  the 

Lagrangian  Function. 

If  we  define  some  vector  5  to  be  an  arbitrary  unit  direction  vector  in  n  space,  then  the  dot 
product  of  Vir(x‘/)  and  S  yields  the  rate  of  change  of  the  objective  function  in  direction  S.  If 
the  dot  product  is  negative  then  the  objective  function  is  getting  smaller  in  that  direction,  and  if 
positive  it  is  increasing.  If  zero,  the  gradient  and  S  are  orthogonal  and  there  is  nothing  to  be  gained 
by  moving  positive  or  negative  in  that  direction.  Applying  the  same  thought  to  the  constraints,  if 
the  dot  product  of  Vzj(x)  and  S  is  zero,  then  S  lies  on  some  iso-potential  surface  of  the  constraint. 
If  the  constraint  is  active  then  this  direction  is  tangent  to  the  constraint  boundary. 

Let  us  assume  a  constraint  is  active  and  direction  S  is  tangent  to  its  boundary.  If  the  dot 
product  of  the  objective  function  and  S  is  nonzero,  the  objective  function  can  be  reduced  by  moving 
along  S.  If  this  is  the  case  then  the  design  cannot  be  a  minimum  for  the  objective  function  can 
be  reduced  without  violating  the  constraint.  However,  if  the  dot  product  is  zero  then  no  move 
from  the  present  design  that  remains  in  the  feasible  region  will  reduce  the  objective  function; 
hence  the  present  design  is  at  least  a  local  minimum.  Also,  since  the  dot  product  of  S  with  the 
gradients  is  zero,  the  gradients  must  be  collinear.  Therefore  an  arbitrary  A  can  be  chosen  such  that 
VF(x")  —  AVz(x1')  =  0,  which  introduces  the  concept  of  the  Lagrangian  function  and  Lagrangian 
multipliers. 

As  Venkayya  (26)  has  demonstrated,  the  constrained  optimization  problem  is  typically  refor¬ 
mulated  with  a  Lagrangian  function  L(x,  A),  defined  as: 

p 

L(x,  A)  =  F(x)  -  A;z;(x)  (8) 

;  =  i 

The  Zj(x)  are  the  p  active  constraints  and  the  A’s  are  the  Lagrangian  multipliers  corresponding  to 
these  constraints.  At  the  minima  the  gradients  of  the  objective  function  and  the  constraints  are 


collinear  and  there  exists  some  collection  of  Lagrangian  multipliers  that  will  nullify  the  sum  of  the 
gradients,  or: 


8L  8F  .  dzj  A  .  , 

dn~  dxi  “  *  -  — 


;= i 


(9) 


The  system  of  Equations  at  9  represent  the  necessary  conditions  of  optimality  and  are  known  as 
the  Kuhn-Tucker  conditions.  If  we  define  a  matrix  e  such  that 


dz, 

e  -  -  — 

_  dF 

then  the  Kuhn-Tucker  conditions  may  be  reformulated  as: 


(10) 


V 

^e0  A_,  =  l  i—  1, . . . ,  n 
;=i 


(11) 


If  the  Kuhn-Tucker  conditions  are  met  and  the  constraints  are  satisfied,  then  the  design  has 
at  least  reached  a  local  minima.  If  the  objective  function  and  inequality  constraints  are  convex, 
and  equality  constraints  are  linear,  then  this  is  also  the  global  minimum  (2:Arora).  Alternatively,  a 
rigorous  search  through  the  feasible  design  space  in  search  of  other  global  minima  candidates  may 
satisfy  the  user  that  the  global  minimum  has  been  reached. 


Methods  For  Moving  About  The  Design  Space 

The  design  optimization  process  requires  moving  through  the  design  space  from  some  x"  in 
search  of  the  x*  that  minimizes  F(x).  Much  work  has  gone  into  developing  algorithms  for  perform¬ 
ing  this  task.  One  common  approach  is  known  as  the  Method  of  Feasible  Directions  (17:Miura). 
Its  robustness  and  efficiency  have  led  to  its  implementation  in  large  structural  analysis  software 
packages  such  as  MSC/NASTRAN  and  ASTROS.  This  method  follows  a  two-phase  approach;  first 
determine  a  search  direction,  then  determine  the  distance  to  step  in  that  direction.  A  less  common 
approach  is  defined  as  the  Optimality  Criteria  Method  (26:Venkayya).  It  is  also  a  two-phase  method 
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which  either  resizes  the  xv  toward  the  optimum  or  scales  it  toward  the  constraint  boundaries.  The 
decision  to  resize  or  scale  is  based  on  the  design  state. 

Method  Of  Feasible  Directions:  In  the  Method  of  Feasible  Directions  (17:Miura),  the  design 
xv  is  iteratively  updated  by  the  expression 

x"+1  =  xv  +  a*S  (12) 

where  5  is  a  unit  pointing  vector  used  to  identify  a  search  direction  and  the  a*  is  a  scaler  move 
parameter  that  defines  the  distance  to  move  along  S.  The  first  critical  step  in  the  optimization  task 
is  to  find  the  search  direction  S.  If  the  design  is  feasible  and  there  are  no  active  constraints,  then 
the  search  direction  that  reduces  the  objective  function  most  rapidly  is  S  =  —  VF.  This  is  known 
as  the  steepest  descent  direction.  Experience  has  shown  that  convergence  is  greatly  enhanced  if 
S  is  modified  towards  the  trend  direction  of  the  last  several  iterations.  If  the  design  is  feasible 
but  there  are  active  constraints,  then  the  steepest  descent  direction  will  likely  violate  a  constraint. 
In  this  case  a  sub-optimization  must  be  performed  to  find  an  S  that  minimizes  VF(*1/)  •  S 
subject  to  the  constraints  Vij(x")  •  5  <  0,  where  the  Zj  are  the  active  constraints  at  the  design. 
If  the  current  design  is  infeasible,  then  a  sub-optimization  must  be  performed  as  above  but  with 
the  objective  function  penalized  and  the  constraint  relaxed  proportional  to  the  magnitude  of  the 
constraint  violation.  This  allows  S  to  point  back  toward  the  feasible  region. 

After  a  suitable  search  direction  has  been  chosen,  the  scalar  step  distance  a*  must  be  de¬ 
termined.  This  is  done  by  a  one  dimensional  search.  By  sampling  the  F(x)  and  zj(  x)  and  their 
gradients  along  S,  interpolation  may  be  used  to  select  an  a*  that  activates  a  new  constraint  or 
minimizes  the  objective  function.  With  a  suitable  step  direction  and  distance  in  hand,  the  design 
is  updated  by  Equation  12. 

At  this  point  the  third  critical  element  of  the  Method  of  Feasible  Directions  comes  into  play; 
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detection  of  convergence  to  the  optimum.  Failure  to  find  a  feasible  search  direction  S  while  in 
possession  of  a  feasible  design  indicates  that  the  Kuhn-Tucker  necessary  conditions  for  optimality 
are  satisfied  and  the  current  design  has  reached  an  optimum.  Another  criteria  is  that  of  having  the 
relative  or  absolute  change  in  the  objective  function  be  below  some  user  specified  tolerance.  The 
optimization  must  also  be  terminated  if  after  a  reasonable  number  of  iterations  a  feasible  design 
cannot  be  found.  If  convergence  is  detected  by  one  of  these  criteria  then  the  optimization  process 
is  stopped.  Otherwise  a  new  search  direction  and  step  parameter  are  sought  and  the  process  goes 


Optimality  Criteria  Method:  In  the  Optimality  Criteria  Method  (2b,  27:Venkayya)  the  design 
xv  is  iteratively  updated  by  one  of  the  following  expressions: 


x"+1  =  x? 


5Z e,]  ^ 

}= 1 


i=  1, 


(13) 


r"+l  =  AjX*'  i-  1, ....  n 


(14) 


Equation  13  resizes  the  sc"  by  stepping  toward  the  optimum  design  and  Equation  14  scales  the 
design  by  stepping  toward  constraint  boundaries.  The  A<  in  Equation  14  is  a  scale  factor  chosen  to 
modify  design  variable  i  so  as  to  bring  the  design  to  a  constraint  boundary. 

The  resizing  scheme  is  based  on  the  Kuhn-Tucker  conditions  expressed  in  Equation  11.  The 
bracketed  sum  at  Equation  13  equals  1  at  the  optimum  design,  indicating  satisfaction  of  the  Kuhn- 
Tucker  conditions  and  producing  no  change  in  xv .  If  the  design  is  not  at  the  optimum  then  the 
bracketed  summation  can  be  thought  of  as  the  ratio  of  the  x‘  at  the  optimum  and  the  current  x". 
Since  the  objective  function  and  constraints  are  typically  nonlinear,  the  relationship  will  not  be 
exact,  but  it  provides  a  factor  that  will  move  the  design  in  the  right  direction.  The  ’’twiddle”  factor 
a  is  used  to  soften  the  magnitude  of  the  move  and  prevent  oscillations.  It  is  adjusted  throughout 
the  optimization  process.  Large  values  of  a  increase  the  number  of  iterations  but  provide  smoother 


10 


convergence.  An  a  <  1  speeds  up  the  iteration  but  may  miss  the  optimum.  Experience  indicates 
that  using  an  a  value  of  2  is  good  in  the  early  iterations  but  to  smooth  convergence  it  should  be 
increased  to  3  or  4  as  the  optimum  is  approached.  For  the  bracketed  summation  to  accurately 
reflect  the  Kuhn-Tucker  conditions  it  would  seem  necessary  to  know  what  constraints  are  active  at 
the  optimum  and  the  values  of  the  Lagrangian  multipliers.  But  Venkayya  (26)  has  noted  that  the 
optimization  process  is  a  series  of  iterations  based  on  approximations,  and  the  approximations  will 
eventually  converge  to  the  appropriate  constraint  set  and  Lagrangian  multipliers.  He  then  details 
a  technique  for  getting  good  first  approximations  for  the  A’s  using  only  gradients  of  the  objective 
function  and  constraints,  weighted  by  each  design  variable’s  contribution  to  the  objective  function. 

The  sealing  scheme  in  Equation  14  moves  the  design  toward  constraint  boundaries  by  applying 
a  scale  factor  A,  to  each  design  variable  Xi,  as  in  Equation  14.  Each  design  variable  may  have  a 
number  of  candidate  scale  factors  to  choose  from  since  each  one  may  in  general  affect  each  constraint, 
and  the  constraints  may  be  violated  or  satisfied  to  differing  degrees.  Venkayya  and  Tischler  (27) 
have  developed  a  compound  scaling  algorithm  that  uses  constraint  function  gradients,  constraint 
function  values,  and  design  variable  values  in  a  weighting  scheme  that  selects  the  proper  factor. 
Two  key  quantities  are  used  in  compound  scaling;  candidate  scale  factors  and  the  relative  sensitivity 
of  the  design  variables  to  the  constraints. 

In  developing  candidate  scale  factors  it  is  important  to  note  that  the  constraints  at  Equations  3 
through  6  are  a  generalization  of  constraints  which  are  typically  expressed  in  some  variation  of: 

Zj(xi,  ...,xn)  <Zj  (15) 

The  zj  is  some  constraint  specified  by  the  user  in  formulating  the  problem,  such  as  a  limitation  on 
stress  or  displacement.  If  we  define  a  parameter  as 

/?>  =  ?-  (16) 
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then  two  candidate  scale  factors  are: 


The  a  is  again  a  "twiddle”  factor  used  to  control  convergence.  Either  scale  factor  may  apply  to  the 
design  variable.  The  sign  of  is  important  in  scale  factor  selection  since  it  indicates  the  direction 
the  constraint  equation  moves  in  response  to  the  variable.  The  sign  of  the  Zj  is  also  important 
for  it  indicates  what  direction  the  constraint  needs  to  move;  whether  from  the  violated  or  feasible 
region.  In  order  to  determine  whether  Ayi  or  Ay  2  is  appropriate  for  a  particular  constraint,  a  new 


parameter  is  defined: 


(19) 


The  sign  of  the  parameter  determines  which  scale  factor  to  use.  If  for  constraint  j  and  design 
variable  i  the  parameter  is  positive,  then  A  j?  is  the  correct  scale  factor,  otherwise  Ayi  should 
be  used.  Applying  this  process  to  all  the  design  variables  leads  to  generation  of  a  scale  factor  table. 
This  is  the  collection  of  candidate  factors  for  each  design  variable  based  on  the  relationship  of  that 
variable  to  each  constraint.  It  now  remains  to  select  the  appropriate  candidate  for  each  design 
variable. 

In  order  to  make  a  proper  selection,  the  design  variables  most  sensitive  to  each  constraint 
must  be  identified.  This  is  done  by  defining  a  new  parameter 


,  _M 

Ui-W 


where  JIJ  is  the  maximum  for  constraint  Zj  .  The  magnitude  of  parameter  mj  provides  a  measure 
of  the  sensitivity  of  constraint  Zj  to  percentage  changes  in  variable  z,.  The  new  parameter  t,} 
normalizes  this  sensitivity  for  each  constraint.  The  design  variables  most  sensitive  to  the  constraints 
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are  identified  by  a  1  and  the  remaining  are  scaled  with  respect  to  these.  The  largest  t,j  associated 
with  design  variable  x,  identifies  the  constraint  Zj  that  is  most  reactive  to  changes  in  the  design 
variable,  and  thus  the  appropriate  scale  factor.  With  this  in  hand,  the  expression  in  Equation  14 
can  be  used  to  update  the  design.  The  entire  scaling  process  is  essentially  bookkeeping.  It  requires 
collecting  the  products  and  ratios  of  gradients,  design  variables,  and  constraint  values  into  tables, 
then  applying  rules  to  those  tables  in  order  to  select  scale  factors. 

The  decision  to  scale  or  resize  is  based  on  the  state  of  the  design.  If  the  design  is  unconstrained 
or  infeasible,  then  it  should  be  scaled  toward  the  constraints.  If  it  is  constrained  then  it  should  be 
resized  toward  the  optimum.  After  updating  the  design  by  Equation  13  or  14,  a  test  for  convergence 
is  applied  as  in  the  Method  of  Feasible  Directions.  This  will  lead  to  termination  of  the  optimization 
process  or  a  new  scale  or  resize  decision. 

Optimization  In  The  Context  Of  Finite  Elements 

The  finite  element  model  of  a  structural  system  is  typically  formulated  as  follows:  Let  x  rep¬ 
resent  the  n  component  vector  of  design  variables  for  the  system.  Let  U  and  P(x )  be  /  component 
vectors  representing  the  generalized  displacements  and  loads  at  nodal  points  of  the  system.  An 
equilibrium  equation  is  formulated  as 


[*■(*)]  U  =  P(x)  (21) 

Where  [K (*)]  is  an  /  by  /  matrix  called  the  stiffness  matrix.  The  stiffness  matrix  and,  in  general, 
the  load  vector  are  functions  of  the  design  variables.  The  displacements  are  obtained  by  .nverting 
the  stiffness  matrix  and  manipulating  Equation  21  to  the  form: 

U  =  [K(x))-1  P(x)  (22) 


13 


Stress,  strain,  and  other  responses  are  then  obtained  by  the  expression 


R  =  [S(z)]  U  =  [S(x)]  [if  (z)]-1  P(x)  (23) 

where  R  is  the  response  vector  of  interest  and  [5(x)]  is  a  response  recovery  matrix  that  relates 
responses  to  displacements.  Like  the  stiffness  matrix,  the  response  recovery  matrix  is  a  property 
of  the  structural  system  and  depends  explicitly  on  the  design  variables,  material  properties  and 
geometry  of  the  system.  Because  of  the  inverted  stiffness  matrix  in  Equations  22  and  23,  the 
explicit  functional  form  of  U  and  R  cannot  generally  be  written;  they  are  implicit  functions  of 
the  design  variables.  Since  constraints  are  placed  on  displacements  and  other  responses,  our  zj  in 
Equations  3  and  4  are  implicit  functions  of  the  design  variables.  This  means  that  to  use  the  finite 
element  model  in  the  optimization  process  requires  the  inversion  of  the  stiffness  matrix  at  each  new 
sample  point,.  This  is  a  very  inefficient  and  impractical  approach  since  the  stiffness  matrix  for  most 
practical  design  problems  is  extremely  large,  and  each  element  may  have  many  design  constraints 
applied  to  it. 

The  optimization  problem  is  made  tractable  by  linearizing  the  finite  element  model  about 
the  current  design  and  ignoring  constraints  that  are  not  likely  to  affect  the  design  at  that  point. 
Efficient  design  sensitivity  analysis  procedures  have  been  developed  for  calculating  derivatives  of 
implicit  functions  with  respect  to  the  design  variables  (2:Arora).  With  values  of  the  objective 
function,  constraints,  and  gradients  at  design  x„,  where  x0  =  {x0\,. .  .  ,x0n),  an  approximate  model 
of  the  optimization  problem  may  be  stated  in  the  form  of  Equations  1  through  6  as  follows: 

Given  the  design  variables 

x  =  (xi,.  (24) 
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and  the  objective  function 


F{x)  =  F(x0)  + 


(xi  -  I„i) 


dFjx) 

dx\ 


"b  *  *  ‘  "b  {Xn  Xon) 


dFjx)' 
dxn  _ 


(25) 


find  the  *  that  will  minimize  the  objective  function,  subject  to  the  following  constraints: 


zj(xo)  + 

zj(x0)  + 


,  dzj(x)  3z;(x)l 

Ui-Xoi)-^ - -  +  •  +(*„  ~x„n 1 

di!  dxn 


,  dzj(x)  3z,(x)l 

(xi-ioi)-^r+"  +(i"'i<’")  dxn 


<0  j  =  l...k 
=  0  j  =  k +  1 . . .  / 


upper 


<0  j  =  1  .  .  .  n 


%j  lower  ^0  i  —  1  .  .  .  M 


(26) 

(27) 

(28) 
(29) 


Since  for  approximate  optimization  we  only  need  the  relative  size  of  the  objective  function,  F(x) 
may  be  simplified  to 


Fr(x)  = 


_  dF( *)  ,  ,  _  dFi *) 

X\— - b - b  In 


dx 


dxn 


(30) 


where  Fr  denotes  the  approximated  relative  objective  function.  The  other  terms  in  Equation  25 
are  constant. 

Using  these  approximations,  optimization  of  structures  modeled  by  finite  elements  is  carried 
out  in  the  following  manner: 


1.  Create  an  approximate  model  at  the  current  design. 

2.  Optimize  the  approximate  model  to  convergence. 

3.  Reevaluate  the  finite  element  model  at  the  new  design. 

4.  Test  for  convergence  at  the  finite  element  level.  Convergence  is  based  on  the  relative  and 
absolute  movement  of  the  objective  function  between  the  finite  element  model  of  step  1  and 
that  of  step  3. 
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5.  If  the  design  has  not  converged  at  the  finite  element  level,  create  another  approximate  model 
and  continue  the  process. 

In  this  way  movement  about  the  design  space  is  driven  by  the  approximate  model.  This  avoids 
many  costly  inversions  of  the  stiffness  matrix  and  renders  the  problem  tractable. 
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III.  ASTROS  AND  MSC/NASTRAN  COMPARED 


MSC/NASTRAN  is  a  large  industry  standard  multidisciplinary  engineering  analysis  tool. 
ASTROS  is  a  relative  newcomer  that  has  challenged  conventional  approaches  to  analysis  by  pio¬ 
neering  the  multidisciplinary  optimization  of  entire  aerospace  systems.  Since  both  are  the  progeny 
of  NASA’s  NASTRAN  analysis  software  package  developed  in  the  1970’s,  they  are  notable  for  their 
similarities.  There  are  also  clear  differences  in  how  each  program  is  directed,  in  how  they  store  and 
manipulate  information,  and  in  their  capabilities. 

Software  Differences 

MSC/NASTRAN  is  described  (15:MacNeal)  as  a  large-scale  digital  computer  program  which 
solves  a  wide  variety  of  engineering  problems  by  the  finite  element  method.  Capabilities  include 
linear  and  nonlinear  static  analysis  with  thermal  loads,  dynamic  structural  analysis,  heat  transfer, 
acoustics,  electromagnetism  and  other  types  of  field  problems.  The  latest  version  includes  a  struc¬ 
tural  optimization  capability  that  can  automatically  modify  the  structural  design  and  associated 
analysis  model  to  satisfy  the  criteria  prescribed  by  the  user  (17:Miura).  Contrary  to  Miura  (17), 
this  capability  was  preceded  by  ASTROS  and  is  therefore  not  unique. 

ASTROS  likewise  is  a  finite  element  based  program  created  to  assist  in  the  design  of  aerospace 
structures  ( 10:Johnson).  Capabilities  include  linear  static  structural  analysis  with  thermal  loads, 
normal  mode  analysis,  aerodynamic  analysis  with  steady  state  aeroelasticity  and  flutter,  and  dy¬ 
namic  analysis  with  transient  and  frequency  response  (18:Neill).  Capabilities  are  clear.y  slanted  to¬ 
ward  the  ’’aero”  disciplines  of  aerospace  and  are  inferior  in  quantity  to  those  available  in  MSC/NASTRAN. 
The  program  performs  multidisciplinary  optimization  within  the  context  of  its  analysis  capabilities. 

The  user  interfaces  with  these  programs  through  an  input  data  file  and  the  subsequent  output 
generated  during  execution.  The  input  data  files  for  MSC/NASTRAN  and  ASTROS  take  a  very 
similar  form.  Both  include  portions  which  control  the  solution  sequence  and  input  the  physical 
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characteristics  of  the  problem  at  hand.  Differences  lie  in  the  logical  arrangement  of  the  solution 
sequence  portion,  the  available  options  for  describing  the  physical  problem,  and  in  the  description 
of  the  optimization  problem. 

Control  of  the  solution  sequence  in  MSC/NASTRAN  is  performed  in  two  sections  of  the  input 
data  file;  the  Executive  section  and  Case  Control  section  (20:Reymond).  The  Executive  section 
identifies  the  job  and  type  of  solution  and  declares  the  general  conditions  under  which  the  solu¬ 
tion  will  be  obtained,  such  as  maximum  time  allowed.  Any  modifications  to  the  norma]  solution 
sequence  are  also  declared  here.  Case  Control  defines  the  subcase  structure  of  the  problem,  such 
as  what  loads  and  constraints  are  to  be  considered.  This  is  where  requests  for  output  are  made.  In 
ASTROS  the  portions  of  the  input  data  stream  that  control  the  solut  ion  sequence  are  the  Executive 
System  Packet  and  Solution  Control  Packet  (18:Neill).  The  input  file  begins  with  a  preface  state¬ 
ment.  to  establish  the  database,  then  the  Executive  System  Packet  calls  a  user-supplied  solution 
sequence  or  modifies  the  standard  one,  if  nece«eai\  li  u&  Executive  is  included  then  the  standard 
sequence  is  executed.  Next,  the  Solution  Control  Packet  defines  the  subcase  structure.  The  most 
notable  cifference  in  this  area  is  in  the  way  uL;i;>hie  „  ,.re  enabled  for  analysis  and  optimization. 
MSC/NASTRAN  has  many  solution  sequences  available  for  the  different  disciplines,  and  these  are 
called  by  declaring  a  solution  number  in  the  Executive.  ASTROS  has  only  one  standard  solution 
sequence  that  encompasses  all  the  available  disciplines.  The  selection  of  disciplines  is  made  by 
simply  rearranging  the  subcase  structure  in  the  Solution  Control  Packet.  MSC/NASTRAN’s  opti¬ 
mization  solution  sequence  SOL  200  is  basically  equivalent,  but  the  ASTROS  approach  is  arguably 
easier  to  understand  and  apply. 

The  physical  problem  being  considered  is  described  within  the  Bulk  Data  section  of  the 
input  data  file.  The  Bulk  Data  entries  available  in  ASTROS  are  essentially  the  same  as  vheir 
counterparts  in  MSC/NASTRAN,  with  the  exception  being  the  description  of  the  design  model  for 
optimization.  MSC/NASTRAN  allows  a  wide  variety  of  property  characteristics  to  be  identified 
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as  design  variables.  ASTROS  allows  only  one  type  of  design  variable  for  each  element  type;  a 
variable  related  to  system  weight,  such  as  thickness  or  area.  MSC/NASTRAN  identifies  a  variety 
of  responses  that  can  be  constraints  or  the  cost  function.  The  responses  that  can  be  used  in 
this  way  are  structural  weight  and  volume;  static  displacement,  stress  and  strain;  internal  force; 
natural  frequency;  buckling  responses  and  composite  lamina  responses  (17:Miura).  ASTROS  has 
Bulk  Data  entries  that  identify  constraints  related  to  aileron  effectiveness,  lift  effectiveness,  flutter, 
modal  frequency,  static  displacement,  stress,  and  strain  (18:Neill).  But  the  program  assumes  all 
optimization  is  based  on  system  weight.  Because  of  this  and  the  limited  design  variables,  the 
Bulk  Data  entries  for  directing  optimization  in  ASTROS  are  simpler  and  incompatible  with  those 
in  MSC/NASTRAN.  The  slight  incren'e  in  complexity  is  a  small  price  to  pay  for  the  improved 
capabilities. 

Other  incompatibilities  in  the  bulk  data  arise  from  MSC/NASTRAN  capabilities  that  AS¬ 
TROS  does  not  support.  Except  for  the  optimization  differences,  an  ASTROS  data  deck  will 
generally  run  in  MSC/NASTRAN.  The  opposite  is  not  true,  however.  MSC/NASTRAN  includes 
numerous  elements,  loads,  parameter  adjusters,  output  requests,  and  other  capabilities  that  are  not 
supported  in  ASTROS.  This  can  limit  the  analyst. 

A  significant  development  in  ASTROS  that  is  invisible  to  the  user  is  the  structure  of  the  data 
base.  The  developers  recognized  a  need  for  a  data  base  that  could  handle  three  distinct  types  of 
data  ( 10:Johnson).  First,  it  had  to  stoie  and  retrieve  large,  sparse  matrices.  The  second  requirement 
is  to  access  individual  data  items  in  entities  such  as  tables,  and  the  third  is  to  access  heterogeneous 
collect  ions  of  unstructured  data  efficiently.  As  no  such  data  base  was  available  commercially,  it  was 
developed  for  ASTROS  and  designated  the  Computer  Automated  Design  Data  Base  or  CADDB. 
The  CADDB  supports  these  different  data  types  and  provides  a  common  structure  for  accessing  all 
three  in  an  efficient  manner.  The  CADDB  also  allows  the  data  base  to  expand  and  contract  with 
the  problem  at  hand  rather  than  resorting  to  fixed  size  arrays. 
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A  notable  similarity  between  MSC/NASTRAN  and  ASTROS  is  their  utilization  of  Vander- 
plaats  MICRO-DOT  algorithm  to  perform  optimization  (17,  10:Miura,  Johnson).  Different  default 
parameters  are  used  by  each  program,  however.  For  instance  MSC/NASTRAN  imposes  a  0.2  rel¬ 
ative  move  limit  while  that  in  ASTROS  is  2.0.  As  a  result,  optimization  in  the  former  may  grind 
along  at  a  much  slower  pace  unless  the  limit  is  relaxed. 

When  both  programs  are  used  side  by  side  on  the  same  problem,  MSC/NASTRAN  is  notable 
for  it’s  robustness  and  flexibility.  Difficulties  that  are  automatically  cared  for  or  worked  around 
in  MSC/NASTRAN  may  crop  up  as  fatal  errors  in  ASTROS.  For  instance,  modal  analysis  of 
flat  plates  in  ASTROS  requires  omission  of  rotational  degrees  of  freedom;  something  that  is  not 
necessary  in  MSC/NASTRAN.  Also,  ASTROS  does  not  support  some  useful  outputs  such  as  nodal 
stress  values,  preventing  the  analyst  from  obtaining  edge  stresses. 

Some  bugs  still  exist  in  the  ASTROS  code.  For  example,  if  the  number  of  constraints  is 
about  ten  times  greater  than  the  number  of  design  variables,  optimization  may  halt  with  a  fatal 
array  dimensioning  error.  When  static  thermal  analysis  is  attempted  on  the  Sun  4  using  either 
Versions  4  or  5,  analysis  terminates  with  a  fatal  error  and  the  system  accuses  the  user  of  defining 
more  than  one  temperature  per  grid  point.  The  same  problems  run  correctly  on  the  VAX.  The 
Flight  Dynamics  Lab  has  confirmed  that  ASTROS  was  poorly  ported  to  the  Sun,  and  is  working  on 
corrections.  Another  interesting  discrepancy  has  to  do  with  Version  5’s  implementation  of  the  four 
noded  quadrilateral  element  QUAD4.  Models  using  this  element  that  run  successfully  in  Version  4 
often  terminate  with  fatal  matrix  singularities  in  Version  5.  The  QUAD4  element  was  modified  for 
Version  5  to  enhance  accuracy  with  composite  materials  (24).  It  appears  that  greater  composite 
material  accuracy  has  driven  the  element  to  singularities  in  other  applications. 
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Comparison  With  Exact  Solutions 

Round  Plate:  Timoshenko  (23)  shows  that  the  center  deflection  w  of  a  uniform  round  plate 
with  clamped  edges  and  loaded  by  a  uniform  pressure  is  given  by 


t'i  = 


qa^_ 

64  D 


(31) 


where  q  is  the  uniform  pressure,  ”a”  the  plate  radius,  and  D  its  flexural  rigidity.  Meirovitch  (16) 
shows  that  the  first  natural  frequency  w  of  such  a  plate  in  Hertz  is  given  by. 


u  =  1.0152 


7T 

2^ 


(32) 


where  p  is  the  area  density  of  the  plate.  The  flexural  rigidity  D  is  given  by 


D  — 


Eh3 

12(1  -  t/2) 


(33) 


where  E  is  Young’s  Modulus,  h  is  the  plate’s  uniform  thickness  and  u  is  Poisson’s  Ratio  for  the 
material. 

A  numerical  experiment  was  performed  letting  "a”  equal  15  inches,  h  equal  0.25  inches,  E 
equal  10700  ksi,  u  equal  0.33,  P  equal  10  psi  and  p  equal  6.475  x  10-5  slug  inches.  Material 
properties  are  nominal  values  for  2024-T3  bare  aluminum  alloy  (5:Denno).  The  theoretical  center 
point  deflection  given  by  Equation  31  is  0.50593  inches  and  the  first  natural  frequency  given  by 
Equation  32  is  111.76  Hertz.  The  same  plate  was  next  evaluated  using  finite  elements.  Figure  1 
illustrates  a  finite  element  model  of  the  plate  and  Figure  2  plots  center  point  deflection  in  inches 
versus  the  number  of  radial  elements  in  the  model.  Both  MSC/NASTRAN  and  ASTROS  converged 
to  within  3.7  percent  of  the  theoretical  result  with  ASTROS  slightly  more  accurate.  Figure  3 
plots  the  first  natural  frequency  of  the  plate  in  Hertz  versus  the  radial  grid  refinement.  Both 
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Figure  1  Model  of  Clamped  Round  Plate  With  Ten  Elements  in  Radial  Direction 
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Figure  2.  Plot  of  Center  Displacement  Versus  Grid  Refinement,  Round  Plate 
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Figure  3.  Plot  of  First  Natural  Frequency  Versus  Grid  Refinement,  Round  Plate 

programs  converged  to  within  0.27  percent  of  the  theoretical  result  with  ASTROS  again  slightly 
more  accurate. 


Plane-strain  Cylinder:  Saada  (21)  shows  that  the  mid-plane  radial  deflection  6  of  a  plane- 
strain  tube  loaded  by  internal  pressure  is  given  by 


Pa7(  1  +  t/)  4  62 

2(6  —  a)E  [(a  +  6)2  "  "  + 


where  P  is  the  internal  pressure,  6  the  outer  radius  and  ”a”  the  inner  radius  of  the  tube.  He  also 
shows  that  the  mid-plane  hoop  stress  an  is  given  by 


a2P  r  462 
°h  ~  62  -  a2  L1  +  (a  +  6)2 


A  numerical  experiment  was  performed  with  "a”  equal  14.875  inches  and  6  equal  15.125 
inches.  This  gave  a  plane-strain  tube  with  a  mid  plane  radius  of  15  inches  and  thickness  of  0.25 
inches.  Using  the  E  and  u  for  aluminum  and  an  internal  pressure  equal  to  10  psi  gives  with 
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Figure  4.  Model  Of  A  Plane  Strain  Tube  With  32  Circumferential  Elements 


Equation  34  a  theoretical  mid-plane  radial  deflection  of  0.00074629  inches.  The  hoop  6tress  is 
594.98  psi  per  Equation  35.  The  tube  was  then  analyzed  with  finite  elements.  Element  boundaries 
were  skewed  with  respect  to  the  ends  of  the  tube  to  convert  any  inter-element  bending  tendencies 
to  membrane  stresses.  Figure  4  illustrates  a  finite  element  model  of  the  tube  and  Figure  5  plots 
mid-plane  displacement  in  inches  versus  the  number  of  circumferential  elements.  MSC/NASTRAN 
converged  to  within  0.2  percent  of  the  theoretical  result  while  ASTROS  only  converged  to  within 
3.7  percent.  Figure  6  is  a  plot  of  the  hoop  stress  o/,  of  the  tube  in  psi  versus  the  number  of 
circumferential  elements.  Both  programs  converge  to  within  about  0.6  percent  of  the  theoretical 
result,  but  MSC/NASTRAN  is  about  0.05  percent  better.  ASTROS  is  slightly  better  on  models 
having  a  coarse  finite  element  grid. 

Clamped  Dome:  In  the  case  of  a  uniform  partial-dome  loaded  by  uniform  pressure  P,  Timo¬ 
shenko  (23)  illustrates  the  general  hypergeometric  series  solution  and  Kraus  (11)  gives  am  example 
of  the  infinite  series  form  of  the  solution  for  the  special  case  of  a  semi-dome.  The  solution  is  an 
infinite  series  of  Gamma  functions  which  becomes  ill-conditioned  and  difficult  to  evaluate  as  the 
thickness  of  the  shell  becomes  small  compared  with  its  radius  (23:Timoshenko).  A  very  good  ap- 
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Figure  5.  Plot  of  Mid-plane  Radial  Displacement  Versus  Grid  Refinement,  Plane-strain  Tube 
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Figure  6.  Plot  of  Hoop  Stress  Versus  Grid  Refinement,  Plane-strain  Tube 
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proximation  can  be  made  by  assuming  a  shape  for  the  shearing  force  that  damps  out  toward  the 
center  of  the  dome.  With  a  radius  to  thickness  ratio  of  30,  Timoshenko  demonstrates  that  the 
approximation  gives  stress  resultants  that  are  good  to  within  10  percent  at  the  extremes  and  follow 
the  theoretical  curves  in  every  respect.  The  approximation  improves  as  the  radius  to  thickness 
ratio  increases.  With  this  approximation  Gibson  (7)  has  solved  for  the  stress  resultant  forces  N<p 
and  Ne  and  the  moments  and  M«,  as  follows: 

N*  —  (36) 

Ne  =  ^(1  -  *V*~O,,l«*(0(*i  -4))  +  sin(/?(a,  -  *))]  -  ~  (37) 

M*  =  — ^-e^»-0l)(sin(/?(oi  -  *))  -  cos(/?(ai  -  *))]  (38) 

Me  =  [sin(^(ai  -  <j>))  -  cos(/?(a,  -  *))]  (39) 

In  these  equations  "a”  is  the  dome  radius  and  m  the  half  angle  subtended  by  the  dome.  The 
parameter  3  is  defined  as: 

/?=</3(TT^)  (40) 

The  0  represents  the  circumferential  direction  and  <f>  the  angle  from  the  dome  center  point,  as  seen 
in  Figure  7  from  Gibson  (7).  The  stresses  within  the  shell  can  be  derived  from  the  stress  resultants 
by  the  following  Equation  (4:Cook): 


<T(.)  = 


yv(.)  12  gM(.) 

h  h3 


(41) 


where  the  *  is  replaced  by  6  or  </>  as  desired,  and  z  is  the  distance  from  the  mid-plane. 

A  numerical  experiment  was  performed  with  "a"  equal  15  inches,  h  equal  0.25  inches,  qj 
equal  35  degrees,  E  and  v  for  aluminum  and  pressure  P  equal  10  psi.  The  radius  to  thickness  ratio 
was  60.  which  is  double  that  used  by  Timoshenko  to  demonstrate  the  quality  of  the  approximate 
solution.  Figure  8  illustrates  a  finite  element  model  of  the  dome.  The  theoretical  values  obtained 
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Figure  8.  Model  Of  A  Clamped  Dome  With  5  Elements  In  The  Radial  Direction 
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ANGLE  PHI  (DEG) 

Figure  9.  Plot  of  Mid-plane  erg  Versus  Angle  4>,  Clamped  Dome 


using  Equations  36  throu'-'  ■  ,vere  compared  with  those  found  by  finite  element  analysis.  Figure  9 
is  a  plot  of  mid-plane  ,r,  ..i  psi  versus  the  angular  distance  from  the  center  point.  Likewise  Figure  10 
plots  (T,p.  The  t'.p  surface  og  is  plotted  in  Figure  11  and  er$  at  the  top  surface  is  shown  in  Figure  12. 
The  stresses  in  each  case  are  in  psi.  MSC/NASTRAN  and  ASTROS  yield  virtually  identical  results 
with  the  former  slightly  better  where  differences  exist.  The  finite  element  solutions  lie  within  about 
9  percent  of  the  theoretical  results  and  their  extreme  value  always  exceeds  theoretical.  This  is  a 
positive  attribute  for  an  analysis  tool  as  the  design  is  driven  by  stresses  not  likely  to  be  exceeded 
in  the  physical  world.  It  is  also  interesting  to  note  that  in  Figure  9  the  finite  element  results  err  in 
the  direction  of  the  exact  hypergeometric  series  solution  illustrated  by  Timoshenko  (23). 

Optimization:  A  numerical  optimization  experiment  was  performed  on  the  64  element  plane- 
strain  tube.  With  the  radial  displacement  from  static  analysis  as  an  upper  constraint,  the  model 
was  optimized  for  minimum  weight  to  see  how  close  MSC/NASTRAN  and  ASTROS  came  to 
the  theoretical  0.25  inch  wall  thickness.  Optimization  started  with  a  thickness  of  10  inches.  With 
default  optimizer  parameters  ASTROS  converged  within  five  iterations,  but  MSC/NASTRAN  took 
fifteen.  This  is  because  the  latter  imposes  a  default  20  percent  move  limit  while  the  former  has  a 
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Figure  11.  Plot  of  Top  Surface  <rg  Versus  Angle  <f>,  Clamped  Dome 
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Figure  12.  Plot  of  Top  Surface  cr^  Versus  Angle  <f>,  Clamped  Dome 


200  percent  limit.  With  the  MSC/NASTRAN  move  limit  relaxed  to  90  percent,  convergence  was 
obtained  in  five  iterations.  Both  programs  converge  to  within  0.02  percent  of  the  theoretical  value. 
The  design  variable  history  is  plotted  in  Figure  13  where  wall  thickness  is  given  in  inches. 
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Figure  13.  Iteration  History  Of  Wall  Thickness,  Plane-strain  Tube 


IV.  OPTIMIZATION  OF  A  DOUBLY  CURVED  SHELL 


The  analysis  of  simple  piate  and  shell  constructions  has  given  confidence  in  the  capabilities 
of  MSC/NASTRAN  and  ASTROS,  and  provided  a  foundation  for  the  analysis  of  more  complex 
structures.  In  this  chapter  optimization  is  probed  more  deeply  in  the  context  of  a  large  nozzle. 

Finite  Element  Model  Of  A  Nozzle 

Physical  Characteristics :  A  finite  element  model  was  created  using  the  deployed  shape  of  the 
Peacekeeper  Stage  ill  nozzle.  1’his  large  nozzle  was  developed  by  the  Missile,  Ordnance  and  Space 
Group  at  Hercules  Corporation.  It  has  a  length  of  78.38  inches,  throat  diameter  of  8.60  inches  and 
exit  diameter  of  69.46  inches  (8:Hercules).  The  nozzle  was  modeled  using  four  noded  quadrilateral 
elements  having  bending,  membrane  and  shear  stiffness.  A  computer  program  was  developed  that 
fits  the  finite  element  mesh  to  the  nozzle  geometry  with  IMSL  cubic  spline  routines  (9),  then 
generates  the  executive,  case  control  and  bulk  data  files  needed  for  analysis  in  MSC/NASTRAN 
and  ASTROS. 

Analysis  time  was  saved  by  exploiting  the  symmetry  of  the  geometry  and  loads.  A  quartered 
model  and  single  element  strip  model  were  used  with  good  results.  A  full  nozzle  with  32  elements 
per  circumference  has  960  elements,  but  only  30  of  these  are  independent.  A  quartered  model 
reduces  the  number  of  elements  to  240,  while  a  strip  model  uses  the  30  independent  elements.  The 
model  generation  program  ha-,  the  flexibility  to  create  either  full,  quartered,  or  strip  models  at  the 
user's  discretion.  It  also  lets  the  user  vary  the  fineness  of  the  mesh.  Figure  14  is  an  illustration  of 
a  quartered  nozzle  and  Figure  15  shows  the  same  nozzle  modeled  by  a  single  strip  of  elements. 

Material  properties  were  those  for  alloy  AISI  687  (13:Lynch);  a  high  temperature,  high 
strength  alloy  of  Nickel,  Chromium,  Cobalt  and  other  trace  elements.  This  alloy  was  chosen  for 
its  yield  strength  at  the  nearly  1700  degrees  fahrenheit  at  the  throat.  AISI  687  has  a  modulus  of 
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Figure  14.  Tinite  Element  Model  Of  Quartered  Nozzle,  Grid  Fineness  Of  32  Elements  Per  Cir¬ 
cumference 


elasticity  of  32400  ksi,  coefficient  of  thermal  expansion  of  7.5X10  6  per  degree  fahrenheit,  and  an 
interpolated  yield  strength  of  68  ksi  at  1700  degrees. 

Pressure  And  Thermal  Loads:  Pressure  loads  were  approximated  using  the  one  dimensional 
Laval  nozzle  equations  for  isentropic  flow.  The  development  of  the  relation  between  pressure  and 
channel  area  is  illustrated  by  Kuethe  and  Chow  (12),  resulting  in  the  equation 


where  the  ”  A"  is  the  nozzle  area  at  the  point  of  interest  and  A’  the  area  at  the  throat.  Pressure 
at  the  point  of  interest  is  denoted  by  P,  and  P0  is  the  chamber  pressure.  The  specific  heat  ratio 
for  the  gas  is  7. 

With  chamber  pressure  and  nozzle  geometry  known,  Equation  42  may  be  solved  numerically 
for  pressure  at  any  axial  location.  The  relation  is  double  valued  except  at  the  throat.  The  correct 
value  is  determined  by  whether  the  point  of  interest  is  up  or  downstream  of  the  throat,  and  the 
pressure  outside  the  nozzle.  For  this  analysis  the  nozzle  was  assumed  to  be  operating  in  a  vacuum, 
thus  the  high  pressure  value  belongs  on  the  chamber  side  of  the  throat  and  the  low  pressure  value 
toward  the  nozzle  exit.  The  model  generation  program  numerically  solves  Equation  42  for  P  using 
the  nozzle  crossectional  area  at  the  element  in  question,  the  throat  area,  a  7  of  1.18  and  a  chamber 
pressure  of  1000  psi.  This  was  applied  to  the  element  as  a  pressure  load.  Element  loads  versus 
location  are  plotted  in  Figure  16. 

The  parameters  used  in  the  development  of  pressure  loads  were  chosen  to  be  as  realistic  as 
possible.  A  1000  psi  chamber  pressure  is  reasonable  for  a  nozzle  of  this  size  (8:Hercules).  The  7  is 
for  a  nitric  acid  oxidizer  and  kerosene  fuel  mixture  used  by  Barrere  (3)  in  an  example  calculation 
of  nozzle  temperatures.  In  a  real  nozzle  the  7  would  not  be  constant,  nor  the  flow  isentropic,  as 
assumed  here.  However,  Equation  42  is  good  first  approximation  to  the  pressure  loads  and  generally 
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Figure  16.  Nozzle  Pressure  Loads 

within  15  percent  of  experimental  values  (6:Galati),  thus  it  provided  a  good  representative  load  for 
nozzle  analysis. 

Temperature  loads  were  modeled  using  the  example  developed  by  Barrere.  His  values  were 
converted  to  degrees  fahrenheit  and  fit  to  this  nozzle’s  geometry.  The  example  was  of  a  nitric  acid 
and  kerosene  fueled  rocket  motor  using  the  nitric  acid  to  regeneratively  cool  the  nozzle.  The  result¬ 
ing  nozzle  wall  temperatures  are  plotted  in  Figure  17.  Although  the  example  is  a  little  outdated,  it 
nonetheless  provided  a  good  representative  temperature  profile  for  analysis  and  material  selection. 


Optimization  Scheme:  Optimization  was  performed  with  the  objective  being  to  minimize 
weight  by  varying  the  wall  thickness.  The  Von  Mises  stress  criteria  was  used  with  a  yield  strength 
of  68  ksi  as  the  constraint.  This  criteria  is  defined  as  (20:Reymond): 


TV  =  \J —  (TzVy  +  rf~  3rIV  (43) 


The  minimum  wall  thickness  was  constrained  to  0.01  inches  and  the  exit  lip  constrained  to  a 
maximum  deflection  of  2  inches  in  any  direction.  Wall  thickness  was  assumed  uniform  in  the 
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Figure  17.  Nozzle  Temperature  Profile 


circumferential  direction,  consistent  with  the  symmetric  loading  and  geometry.  ASTROS  Versions 
Four  and  Five  do  not  support  nodal  temperature  loads  on  the  SUN  4,  therefore  it  was  run  without 
thermal  loads.  MSC/NASTRAN  was  run  with  and  without  thermal  loads  for  comparison  purposes. 

The  Modified  Feasible  Directions  search  scheme  was  used  by  MSC/NASTRAN  and  AS¬ 
TROS.  Both  programs  employ  the  MICRO-DOT  algorithm  for  optimization  but  initiate  it  with 
different  default  optimizer  parameters.  Of  the  45  parameters  in  the  algorithm  (25:Vanderplaats), 
MSC/NASTRAN  provides  for  the  easy  modification  of  13  and  ASTROS  documents  26.  Both  pro¬ 
grams  allow  the  experienced  user  to  vary  all  the  parameters.  Numerous  optimization  runs  were 
performed  to  compare  the  goodness  of  the  default  parameters  and  identify  those  which  have  the 
best  effect  on  the  optimization  process. 


Results 

The  strip  model  ran  well  in  MSC/NASTRAN  with  adjustment  of  the  optimizer  parameters, 
but  did  not  converge  with  defaults.  Figure  18  compares  the  optimization  progress  with  different 
parameters.  The  dramatic  improvement  in  the  first  step  was  obtained  by  changing  the  move  limit 
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Figure  18.  Comparison  Of  Optimization  Progress  Using  Default  And  Modified  Optimization  Pa¬ 
rameters  In  MSC/NASTRAN,  Nozzle  Model  With  Grid  Fineness  Of  32  Elements  Per 
Circumference 


parameter  DELP  from  its  default  of  0.2  to  2.0.  This  resulted  in  a  slight  overshoot  that  settled  in 
four  iterations  and  converged  within  ten.  Adjusting  the  convergence  detection  parameters  CONV1 
and  CON  V2  provided  a  further  improvement  in  the  objective  of  about  2  percent.  These  parameters 
determine  the  relative  and  absolute  change  in  the  objective  for  convergence.  The  default  values  are 
0.001  and  0.01  respectively.  They  were  both  changed  to  0.1X10-12 

The  strip  model  and  the  quartered  mode!  were  used  in  ASTROS.  Figure  19  compares  the 
optimization  progress  between  the  two  models.  Although  the  ASTROS  optimizer  default  values 
were  an  improvement  over  those  in  MSC/NASTRAN,  the  design  never  reached  the  optimum. 
Varying  the  optimizer  parameters  had  little  effect.  Convergence  detection  parameters  were  adjusted 
to  0.1X10“ 12  and  0.1X  10" 18,  and  the  design  variabes  permitted  to  step  to  one  thousandth  of  their 
initial  value  on  the  first  step.  But  the  design  still  changed  at  the  same  rate  and  ground  to  a  halt 
short  of  the  optimum.  The  quartered  model  was  then  run  and  it  converged  at  the  same  rate  but 
settled  to  the  MSC/NASTRAN  optimum  in  eleven  iterations.  Apparently  ASTROS  has  difficulty 
optimizing  extremely  light  structures  without  modification  of  some  yet  undiscovered  parameter. 
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Figure  19.  Comparison  Of  Optimization  Progress  Using  Default  And  Modified  Optimization  Pa¬ 
rameters  In  ASTROS,  Nozzle  Model  With  Grid  Fineness  Of  32  Elements  Per  Circum¬ 
ference 

Moving  to  an  equivalent  but  more  massive  model  overcame  this  difficulty. 

The  best  MSC/NASTRAN  and  ASTROS  optimization  runs  are  compared  in  Figure  20.  Both 
programs  converged  to  final  values  that  agree  to  within  1.3  percent.  However,  MSC/NASTRAN 
arrived  in  the  neighborhood  of  the  optimum  much  more  quickly  than  ASTROS.  The  resulting 
thickness  distribution  is  plotted  in  Figure  21  with  thickness  is  in  inches.  The  programs  are  in  good 
general  agreement  and  reach  the  0.01  inch  minimum  thickness  constraint  at  about  20  inches  axially. 
The  distributions  at  the  throat  reflect  the  rapidly  changing  loads  and  geometry.  MSC/NASTRAN 
has  more  sharply  defined  thickness  variations  in  the  throat  area  resulting  from  the  fourteen  iter¬ 
ations  allowed  it  near  the  optimum.  In  designing  a  real  nozzle,  the  peak  values  from  an  analysis 
such  as  this  would  be  used  to  establish  a  smoothly  varying  thickness  distribution  from  chamber  to 
exit,  with  some  margin  of  safety  added.  Despite  the  discrepancies  at  the  throat,  both  programs 
would  provide  good  data  for  such  a  distribution.  If  ASTROS  were  used,  however,  there  would  be 
some  risk  of  exceeding  margins,  since  the  peak  values  in  ASTROS  are  consistently  below  those  in 
MSC/NASTRAN 
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Figure  20.  Comparison  Of  Optimization  Progress  Using  MSC/NASTRAN  and  ASTROS,  Nozzle 
Model  With  Grid  Fineness  Of  32  Elements  Per  Circumference 


Figure  21  Comparison  Of  Final  Thickness  Distribution  Using  MSC/NASTRAN  And  ASTROS, 
Nozzle  Model  With  Grid  Fineness  Of  32  Elements  Per  Circumference 
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Figure  22.  Comparison  Of  Major  Principle  Stress  Distribution  Using  MSC/NASTRAN  and  AS¬ 
TROS,  Nozzle  Model  With  Grid  Fineness  Of  32  Elements  Per  Circumference 

The  resulting  major  and  minor  principle  stresses  are  plotted  in  Figures  22  and  23.  Once 
again  both  programs  are  in  very  good  agreement  with  the  exception  of  an  interesting  deviation  at 
the  throat.  Major  principle  stresses  are  in  the  circumferential  direction  and  minor  stresses  are  in 
the  axial  direction.  The  circumferential  stress  is  generally  in  tension  since  the  radius  of  curvature 
lies  on  the  pressure  side  of  the  structure.  However,  at  the  throat  the  axial  radius  of  curvature 
lies  outside  the  nozzle.  This  causes  a  large  compressive  spike  in  the  axial  distribution  which  the 
geometry  couples  into  the  circumferential  direction  as  a  trend  toward  compression.  This  is  seen 
in  both  the  MSC/NASTRAN  and  ASTROS  stress  distributions.  The  ASTROS  model  is  about 
0.01  inches  thicker  in  the  vicinity  of  the  throat,  and  this  support  keeps  the  element  from  going 
into  outright  compression.  The  MSC/NASTRAN  model  is  thinner,  so  the  throat  does  go  into 
compression.  Despite  these  variations,  the  Von  Mises  stress  criteria  are  satisfied  in  both  cases, 
and  the  stresses  have  been  shown  to  be  reliable.  The  discrepancy  at  the  throat  is  a  result  of 
slightly  different  thickness  distributions,  and  these  in  turn  are  a  product  of  the  differences  in  the 
optimizer  parameters.  The  nozzle  was  also  optimized  using  a  finite  element  model  with  twice  the 
grid  refinement  .  There  was  little  change  in  the  thickness  and  stress  distributions. 
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Figure  23.  Comparison  Of  Minor  Principle  Stress  Distribution  Using  MSC/NASTRAN  and  AS¬ 
TROS,  Nozzle  Model  With  Grid  Fineness  Of  32  Elements  Per  Circumference 


The  nozzle  was  optimized  in  MSC/NASTRAN  with  static  thermal  loads.  This  was  done  to 
quantify  the  thermal  stresses  being  overlooked  because  of  ASTROS’  inability  to  support  thermal 
analysis  on  the  Sun  4.  Stresses  were  low  and  the  optimization  quickly  moved  the  thickness  dis¬ 
tribution  to  the  0.01  minimum  throughout.  The  resulting  major  principle  stress  distribution  is 
plotted  at  Figure  23.  With  the  exception  of  a  3200  psi  tensile  load  at  the  throat,  the  loads  are 
minimal  and  of  a  few  hundred  psi.  The  throat  load  is  the  result  of  the  rapid  change  in  nozzle 
diameter  and,  interestingly,  opposes  the  stress  resulting  from  the  pressure  load.  When  the  model 
was  optimized  in  MSC/NASTRAN  with  both  pressure  and  thermal  loads,  the  final  objective  was 
about  0.2  percent  smaller  than  that  obtained  with  pressure  loads  alone.  This  reflects  the  balancing 
of  opposing  stresses  at  the  throat,  and  could  lead  to  a  thickness  distribution  that  is  inadequate  at 
engine  ignition  when  the  nozzle  is  cool.  In  this  application  the  inclusion  of  thermal  loads  does  not 
significantly  improve  the  analysis.  The  real  utility  in  knowing  the  thermal  distribution  is  in  the 
material  selection  process,  where  use  of  high  temperature  yield  values  is  critical.  Thermal  analysis 
should  be  done,  and  it  is  unfortunate  that  ASTROS  was  unsupportive,  but  good  results  can  be 
obtained  nonetheless. 
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Figure  24. 


Principle  Stress  Distribution  Resulting  From  Thermal  Loads  Using  MSC/NASTRAN, 
Nozzle  Model  With  Grid  Fineness  Of  32  Elements  Per  Circumference 
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V.  AXISYMMETRIC  FINITE  ELEMENT  ANALYSIS  IN  ASTROS 


Axisymmetric  shell  finite  elements  offer  an  alternative  to  the  analysis  of  axially  symmetric 
shells  with  flat  elements.  We  have  seen  the  results  that  may  be  obtained  with  flat,  quadrilateral 
elements.  This  chapter  will  discuss  the  development  of  a  Mindlin  type  axisymmetric  shell  element 
and  illustrate  its  use  in  ASTROS. 


Element  Development 

The  following  development  carries  through  the  steps  alluded  to  but  not  explicitly  demon¬ 
strated  by  Cook  (4)  in  his  development  of  a  Mindlin  axisymmetric  shell  element.  Element  geometry 
and  coordinates  are  illustrated  in  Figure  25.  As  in  the  figure,  t  denotes  element  thickness  and  z 
is  the  distance  from  the  midplane.  The  length  of  the  element  is  L  and  the  degrees  of  freedom 
supported  at  each  node  ring  are  u,  w  and  0.  The  u  and  w  may  be  rotated  to  the  global  coordinates 
U  and  W.  The  angle  the  element  makes  with  global  coordinate  U  is  denoted  4>.  The  circumferential 
direction  is  not  shown  but  is  denoted  by  0. 

The  element  nodal  displacement  vector  de  is  configured  as: 


dc  = 


ivi  v  1 


01  U’2  U2  02 


T 


(44) 


In  the  axisymmetric  case  these  displacements  refer  to  nodal  rings  rather  than  the  more  familiar 
nodal  points.  Because  of  this,  elements  containing  nodal  points  cannot  be  attached  to  axisymmetric 
elements. 

Using  the  generalized  coordinate  £  where  £  =  the  shape  functions  are: 


Ni=\(l~0  (45) 

N2  =  \(l+0  (46) 
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0  Ni  0  0  N2  o' 

[TV]  =  Nx  0  N3  N2  0  -N3  (48) 

0  0  JV[  0  0  N2 

If  we  define  an  operator  matrix  [D]  as 


(49) 


then  the  element  strains  are  obtained  by  the  operation: 


c  =  [D][N]d  =  [B)d  (50) 

This  gives  e  in  the  form  e  =  f t  (g  Kf  Kf  ytt  J  where  the  e  and  7  are  strains  and  the  k 

are  curvature  changes. 
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The  material  property  matrix  is  denoted  by  [E]  and  defined  as: 


Et  uEt 
uEt  Et 


[E}=  0  0 

0  0 

0  0 


0  0 

0  0 

E  t3  vEt* 


fEt1 
12(1  — i/2) 

0 


0 

5 Et 

12(1  +  0 


Premultiplying  e  by  [ E ]  returns  membrane  forces  A/(»),  and  bending  moments  M(.\  in  the 
shell.  The  (*)  is  replaced  by  s  or  6  depending  on  which  stress  resultant  is  desired.  Using  the 


A^.)  12zA/(.) 


stresses  are  obtained  through  matrix  [S>]  of  the  form: 


'  7  0  U  00" 

0  \  0  £  0 

±  0  0  0  0 

[S,]  =  0  ±  0  0  0  (53) 

±  0  0  0 

0  ±  0  ^  0 

0  0  0  o  1  _ 

If  postmultiplied  by  the  element  strain  vector,  this  matrix  will  return  stresses  in  the  element 
coordinate  system.  These  stresses  are  at  the  outer  surface,  midplane  and  inner  surface  of  the  shell. 
The  form  of  the  stress  vector  returned  is: 


J  mid  pta  ne 

^.„«r 

. 


With  these  matrices  the  element  stiffness  matrix  [fc]  is  formed  by  the  following  integration: 


[*]  =  \  /’  [B]T[^][B]^ 
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The  matrix  is  then  rotated  to  the  global  coordinate  system  for  assembly  into  a  global  stiffness 
matrix. 

Loads  in  general  do  not  fall  at  node  rings.  A  vector  of  equivalent  nodal  loads  /tn;  is  formed 
for  pressure  load  P(£)  on  the  element  by  the  integration: 

/e»i  =  |  J1  [N}TP(Odt  (56) 

After  the  finite  element  solution  has  been  obtained,  element  stresses  are  extracted  from  the 
global  displacement  vector  dg  as  follows: 

<r-  \S9]d,  =  [Se][E}[B}{R)dg  (57) 

The  [J?j  is  a  transformation  matrix  that  rotates  the  global  displacements  dg  to  displacements  in  the 
element  coordinate  system.  [5S]  is  the  stress  recovery  matrix  which  converts  global  displacements 
to  element  level  stresses. 

ASTROS  Implementation 

The  Mindlin  \xisymmetric  element  was  implemented  in  ASTROS  by  externally  programming 
the  necessary  matrices,  creating  a  pseudo  axisymmetric  problem  in  ASTROS  with  rod  elements, 
then  overwriting  the  internal  matrices  with  those  created  externally.  The  program  discussed  in 
Appendix  R  was  developed  to  create  the  global  stiffness  matrix,  load  vector,  and  stress  recovery 
matrix.  The  integrations  required  by  Equations  55  and  56  are  performed  using  a  seventh  order 
Gaussian  Quadrature  scheme  with  weights  and  function  evaluation  points  provided  by  Stroud  and 
Secrest  (22)  The  element  level  matrices  thus  obtained  are  rotated  to  the  global  coordinate  system 
and  assembled  into  global  level  matrices.  Likewise  the  element  level  stress  recovery  matrix  is 
obtained  by  performing  the  multiplications  of  Equation  57.  These  are  assembled  into  a  global  level 
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stress  recovery  matrix.  The  matrices  are  then  put  in  a  direct  matrix  input  form  for  bulk  data  er  ry 
into  ASTROS. 


To  flush  out  the  unwanted  matrices  and  load  in  the  externally  created  ones,  ASTROS  requires 
modification  of  the  standard  MAPOL  solution  sequence  at  the  executive  level.  For  ASTROS  version 
5  this  was  done  through  the  following  edit  sequence: 


INSERT  88 

MATRIX  [MYKGG] , [MYPL] , [MYSREC] , [STROUT]  ; 
INSERT  1613 
[KGG]  :=  [MYKGG]; 

INSERT  1638 
[PG]  :  =  [MYPL]  ; 

REPLACE  2303,2311 
[STROUT]  : =  [MYSREC]  *  [UG] ; 

CALL  UTMPRT(0 , [STROUT] , [UG] ) ; 


The  insertion  at  line  88  creates  matrix  entities  MYKGG,  MYPL,  and  MYSREC  within  AS¬ 
TROS.  These  are  the  global  stiffness,  load,  and  stress  recovery  matrices  to  be  loaded  via  direct 
matrix  input.  The  STROUT  matrix  outputs  stresses.  At  the  appropriate  point  the  ASTROS  global 
stiffness  matrix  KGG  and  load  vector  PG  are  replaced  by  the  externally  generated  substitutes.  So¬ 
lution  occurs  normally,  albiet  with  radically  modified  matrices.  The  normal  output  sequence  from 
lines  2303  through  2311  is  no  longer  useful,  and  is  replaced  with  a  matrix  multiplication  to  produce 
stresses  and  an  output  command  to  directly  print  stresses  and  displacements. 

All  that  remained  was  to  create  a  pseudo  axisymmetric  bulk  data  problem  with  CONROD 
elements.  This  arranged  the  solution  sequence,  constraints,  matrix  sizes  and  the  internal  structure 
of  ASTROS  for  solution  of  the  externally  generated  matrices.  CONROD  is  a  rod  finite  element 
that  has  an  element  stiffness  matrix  and  degrees  of  freedom  of  similar  size  and  compatible  with 
the  axisymmetric  finite  element.  Because  of  this  it  is  useful  for  creating  the  pseudo  problem.  The 
grid  points  must  match  the  location  of  the  nodal  rings  used  to  create  the  stiffness,  load,  and  stress 
recovery  matrices,  and  they  must  lie  in  the  X-Z  plane.  The  Z  coordinate  in  ASTROS  matches 
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the  global  U  coordinate  a1  _.ng  the  axial  direction  in  the  axisymmetric  problem.  Likewise  the  X 
coordinate  matches  the  global  W.  Constraints  and  CONROD  elements  are  applied  to  these  grid 
points  as  they  would  be  in  the  axisymmetric  problem.  Loads  and  material  property  entries  are 
necessary  for  proper  internal  sizing,  but  can  contain  dummy  values  since  externally  generated  data 
will  take  over  before  the  solution  begins. 

Results 

A  numerical  experiment  was  performed  using  the  geometry  and  properties  of  the  plane  strain 
tube  used  in  Chapter  3.  The  direct  matrix  input  data  created  for  the  axisymmetric  problem  were 
merged  with  a  CONROD  pseudo  model  of  the  tube  and  the  appropriate  executive  modifications. 
The  theoretical  mid-plane  radial  deflection  is  0.00074629  inches  and  theoretical  hoop  stress  is  594.98 
psi.  The  ASTROS  axisymmetric  solution  returned  a  radial  deflection  of  0.00074952  inches  and  hoop 
stress  of  600.00  psi.  This  is  an  error  of  0.43  percent  for  displacement  and  0.84  percent  for  stress, 
comparing  favorably  with  the  0.6  percent  obtained  using  32  quadrilateral  flat  elements. 

This  problem  was  formulated  with  a  12  by  12  matrix  containing  the  6  by  6  axisymmetric 
stiffness  matrix.  In  contrast,  the  quadrilateral  element  problem  required  solution  of  a  1932  by 
’932  stiffness  matrix  for  the  full  tube  or  a  24  by  24  matrix  if  symmetry  was  fully  exploited.  The 
axisymmetric  element  offered  comparable  accuracy  while  reducing  the  size  of  the  problem,  but  in 
this  formulation  disallowed  consideration  of  non-axisymmetric  loads. 

The  optimized  32  element  per  circumference  nozzle  model  of  Chapter  4  was  also  reformulated 
as  an  axisymmetric  problem.  This  produced  an  extremely  large  bulk  data  deck  which  ASTROS  was 
not  successful  in  solving.  In  describing  the  bulk  data  entry  for  direct  matrix  input,  the  ASTROS 
User’s  Manual  ( 1 8: Neill )  implies  that  only  nonzero  entries  need  be  loaded  In  fact,  if  the  entire 
matrix  is  not  loaded  including  all  the  zeros,  a  number  of  strange  errors  crop  up  in  the  solution 
sequence.  In  both  versions  4  and  5  the  program  typically  accuses  the  user  of  trying  to  load  beyond 
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the  dimensioned  size  of  the  matrix,  or  of  improperly  using  complex  arguments.  These  problems 
disappeared  when  the  entire  matrix  was  loaded,  but  this  drove  the  size  of  the  direct  matrix  entry 
to  the  order  of  10,000  lines.  This  resulted  in  input/output  problems  on  the  Sun  4. 
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VI.  CONCLUSIONS 


The  Utility  Of  MSC/NASTRAN  And  ASTROS 

MSC/NASTRAN  and  ASTROS  provide  accurate  analysis  and  optimization  results  on  plate 
and  shell  structures,  providing  the  finite  element  model  is  carefully  constructed  with  a  fine  grid 
structure.  Results  on  simple  shells  have  been  demonstrated  to  within  10  percent  of  theoretical 
values,  and  responses  obtained  for  a  large  nozzle  make  good  physical  sense.  From  an  analysis 
standpoint  the  differences  between  the  programs  are  negligible.  Optimization  results  are  similar, 
but  each  program  steps  down  differing  iteration  paths  driven  by  different  optimizer  parameters. 

MSC/NASTRAN  is  more  robust  with  a  large  element  set  and  many  other  options  available 
to  the  analyst.  It  has  well  developed  documentation  suitable  for  both  new  and  experienced  users. 
Optimizer  parameters  are  clearly  documented  and  can  greatly  speed  the  optimization  process.  The 
program  is  capable  of  analyzing  thermal  loads  on  the  Sun  4.  Software  bugs  were  not  encountered. 

ASTROS  is  hindered  by  a  few  software  errors.  If  the  ratio  of  constraints  to  design  variables 
exceeds  about  ten,  the  optimizer  is  driven  into  an  array  dimensioning  error.  The  program  as 
ported  to  the  Sun  Workstation  has  difficulties  with  thermal  analysis.  The  documentation  available 
f'  •  ASTROS  is  also  a  limiting  factor.  As  a  whole  the  documentation  is  accurate  but  not  particularly 
user  friendly.  The  standard  complement  of  manuals  is  available,  but  these  do  not  accurately  reflect 
the  newer  versions  of  ASTROS.  System  Release  .  ates  must  be  used  in  conjunction  with  the  manuals 
to  obtain  correct  information  in  some  cases.  The  small  element  set  and  limited  capabilities  can 
also  hinder  the  analyst. 

On  the  positive  side,  with  the  capabilities  it  has  ASTROS  offers  cost  free  equivalence  to 
MSC/NASTRAN.  Software  bugs  are  being  actively  pursued  and  corrected,  and  an  optimality  cri¬ 
teria  optimization  scheme  is  being  implemented.  The  inclusion  of  this  optimality  criteria  seaml) 
technique  in  the  next  release  should  significantly  reduce  the  consumption  of  computer  time  on  large 
scale  problems. 
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The  Optimization  Process 


Both  MSC/NASTRAN  and  ASTROS  use  the  Method  of  Feasible  Directions  for  optimization. 
This  is  a  robust  technique  that  moves  about  the  design  space  based  upon  gradient  information. 
Movement  is  meticulous  and  deliberate.  A  sub-optimization  must  be  performed  in  order  to  deter¬ 
mine  the  direction  to  move.  A  one-dimensional  line  search  must  be  performed  to  determine  the 
magnitude  of  the  movement.  This  requires  sampling  the  design  space  and  constraints  at  several 
points  and  interpolating  based  on  the  results.  The  method  is  well  conditioned  but  its  utility  is 
limited  to  about  300  design  variables,  sis  are  all  other  popular  methods. 

The  Optimality  Criteria  Method  also  moves  in  the  design  space  based  on  gradient  information, 
but  gradients  are  used  differently.  The  design  leaps  about  the  design  space  driven  by  simple  ratios 
of  the  design  information,  modified  by  ’’twiddle”  factors.  No  sub-optimizations  or  one  dimensional 
searches  are  required.  The  method  is  sensitive  to  ’’twiddle”  factors  and  may  be  unstable  if  they 
are  not  just  right.  On  the  other  hand,  moving  about  the  design  space  is  less  computationally 
burdensome  and  less  sensitive  to  the  number  of  design  variables  in  the  problem.  This  method 
should  move  optimization  beyond  present  design  variable  limits. 

Nozzle  Analysis  And  Optimization 

The  doubly  curved  nozzle  shape  can  be  adequately  analyzed  in  both  MSC/NASTRAN  and 
ASTROS.  Key  areas  of  consideration  are  thermal  effects  and  stress  distributions.  The  analyst 
may  also  be  concerned  with  displacements,  but  these  arc  likely  to  be  of  little  consequence.  The 
displacements  seen  in  the  analysis  done  here  were  of  the  order  of  a  tenth  of  an  inch  at  the  most. 

Nozzle  temperatures  are  extremely  high  and  constrained  by  material  limitations.  Because 
of  the  thermodynamics  of  nozzle  flow,  the  highest  temperatures  are  near  the  throat.  Material 
properties  may  be  strongly  temperature  dependent,  particularly  yield  strength.  For  instance,  the 
room  temperature  yield  strength  of  alloy  AISI  687  (13: Lynch )  is  double  that  at  1700  degrees 
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fahrenheit.  If  the  room  temperature  yield  strength  were  used  in  optimization  at  the  1700  degree 
throat,  catastrophic  failure  would  result.  Therefore  high  temperature  material  properties  must  be 
used  for  nozzle  analysis. 

Static  thermal  loads  also  induce  stresses,  but  these  have  been  shown  to  be  lower  order  effects 
for  the  nozzle  considered  here.  Thermal  loads  using  composite  materials  and  transient  or  erosive 
effects  could  be  significant,  but  are  beyond  the  scope  of  this  discussion.  Thermal  stresses  arise  in 
shells  from  nodal  constraints,  changes  of  curvature,  and  thermal  gradients.  All  three  of  these  effects 
are  strong  in  the  vicinity  of  the  throat  and  can  affect  optimum  thickness  distributions.  Thermal 
stresses  at  the  throat  oppose  those  which  are  pressure  induced,  leading  to  a  thinner  optimum  wall. 
The  thickness  may  be  adequate  at  operating  temperatures  but  insufficient  for  loads  at  ignition 
when  the  nozzle  is  cold.  Nozzles  should  be  analyzed  with  and  without  thermal  effects  to  guard 
against  this. 

Pressure  induced  stresses  are  more  critical  than  those  induced  thermally.  Four  distinct  stress 
distribution  regions  can  be  identified  and  linked  to  the  geometry  and  loading  of  the  nozzle.  Region 
I  is  where  the  nozzle  terminates  at  the  chamber.  Region  II  is  the  throat  area.  Region  III  lies  in 
the  high  pressure  region  aft  of  the  throat.  Region  IV  is  the  lower  pressure  region  aft  of  the  throat 
where  stress  constraint*  are  no  longer  active. 

Because  of  the  nature  of  supersonic  nozzle  flow,  pressures  rapidly  drop  aft  of  the  throat.  Due 
to  this  drop  in  pressure,  the  nozzle  wall  thickness  toward  the  exit  is  no  longer  constrained  by  stress 
but  is  driven  by  other  considerations,  such  as  manufacturability,  buckling,  and  dynamic  loads.  For 
the  nozzle  evaluated  here,  the  optimization  process  quickly  moved  wall  thickness  aft  of  20  inches 
to  the  hypothetical  manufacturability  constraint  of  0.01  inches.  This  area  is  denoted  as  region  IV. 
The  dominant  radius  of  curvature  is  in  the  circumferential  plane  on  the  pressure  side,  making  the 
dominant  stress  a  tensile  hoop  «tress  in  the  circumferential  direction.  The  wall  is  tilted  with  respect 
to  the  axial  direction,  so  there  is  a  cumulative  compressive  axial  direction  stress.  But  pressures  are 
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low  and  the  angle  of  the  wall  with  the  axial  direction  is  small,  making  this  a  secondary  effect. 


Stress  distributions  near  the  throat  are  much  more  interesting  and  unpredictable.  As  the 
throat  is  approached  from  the  aft  end,  pressures  rise  and  eventually  deactivate  the  side  constraints, 
signaling  entry  into  region  III.  Wall  thickness  is  driven  by  a  stress  constraint  such  as  the  Von  Mises 
stress  criteria.  Hoop  stress  varies  proportional,  and  axial  stress  inversely  proportional  to  radius. 
Therefore  in  region  III  a  tradeoff  occurrs  between  axial  and  hoop  stress  to  keep  the  Von  Mises 
stress  maximized.  This  is  seen  in  Figures  23  and  24  where  the  circumferential  or  major  principle 
stress  decreases  toward  the  throat,  while  the  axial  or  minor  principle  stress  increases  in  magnitude. 

Continuing  up  the  nozzle,  the  shell  curvature  inflects  to  form  the  throat,  signaling  the  onset 
of  region  II.  In  this  area  the  radius  of  curvature  in  the  axial  plane  approaches  and  may  exceed  that 
in  the  circumferential  plane.  The  change  in  curvature  superimposes  a  compressive  hoop  stress  on 
the  large  compressive  axial  stress.  The  total  axial  stress  dominates  and,  because  of  the  geometry, 
imparts  a  compressive  circumferential  stress.  This,  combined  with  the  circumferential  hoop  stress 
reduction  due  to  decreasing  diameter,  leads  to  the  trend  reversal  seen  in  circumferential  stress  at 
the  throat  .  The  reversal  allows  a  thinning  of  the  nozzle  wall  in  this  region. 

At  region  I  there  is  a  dramatic  axial  stress  reduction  due  to  the  high  pressure  load  component 
opposing  the  cumulative  effects  from  the  rest  of  the  nozzle.  Radial  stresses  again  dominate,  and 
there  is  a  considerable  thickening  of  the  nozzle  wall  to  support  the  high  pressures. 

Artsymmelnc  Elements 

Axially  symmetric  problems  such  as  the  nozzle  can  be  modeled  using  axisymmetric  finite 
elements.  The  Mindlin  axisymmetric  element  is  nearly  as  accurate  as  a  four  noded  quadrilateral 
element  model  with  a  grid  fineness  of  160  elements  per  circumference.  The  axisymmetric  model 
tested  here  uses  a  6  by  6  stiffness  matrix  to  solve  the  problem  posed  by  a  1932  by  1932  quadrilateral 
element  stiffness  matrix  This  is  a  significant  simplification.  The  matrix  flushing  scheme  provided 
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a  quick  test  of  the  Mindlin  element’s  accuracy  and  demonstrated  ASTROS’  flexibility,  ease  of 
modification,  and  potential  for  the  implementation  of  such  an  element. 

The  axisymmetric  element  considered  here  is  not  attachable  to  elements  with  nodal  points 
and  does  not  support  non-axisymmetric  loads.  The  loading  difficulty  can  be  overcome  by  Fourier 
analysis  (14:MacNeal).  The  inability  to  attach  to  other  element  types  is  a  disadvantage  common 
to  axisymmetric  elements,  precluding  their  use  in  large  complex  models. 

Lessons  Learned  In  The  Optimization  Of  Shells 

The  following  general  guidelines  have  been  gleaned  from  this  investigation  and  will  prove 
useful  to  the  analyst  seeking  the  "correct”  answer  to  a  shell  problem: 

1.  The  analyst  should  not  be  deceived  by  the  remarkably  accurate  results  finite  element  analysis 
gives  for  flat  plates  and  singly  curved  shells  under  simple  loading.  A  complex,  doubly  curved 
shell  is  not  likely  to  give  results  with  an  error  less  than  the  ±  10  percent  seen  for  the  clamped 
dome. 

2.  Thermal  loads  will  induce  stresses  in  response  to  constraints,  curvature  changes,  and  thermal 
gradients.  Material  properties  used  in  analysis  must  also  take  temperatures  into  account. 

3.  Optimization  should  be  driven  to  convergence.  Unless  the  plot  of  the  objective  or  cost  function 
asymptotically  approaches  some  value,  the  optimum  has  not  been  reached.  This  is  seen 
in  Figure  19  where  ASTROS  terminated  optimization  before  reaching  the  minimum  wall 
thickness,  but  claimed  convergence. 

4.  Optimizer  parameters  should  be  varied.  Accepting  default  values  may  give  very  slow  conver¬ 


gence. 


5.  Grid  fineness  should  be  varied.  Accurate  results  are  obtained  when  the  results  are  no  longer 
affected  by  increasing  grid  refinement.  This  was  seen  in  Figures  2,  3,  5  and  6  with  the  round 
plate  and  tube. 

6.  Formulating  the  problem  with  a  different  finite  element  model  may  be  helpful  and  collaborate 
results.  Going  to  a  quartered  model  was  the  only  way  to  get  good  optimization  results  in 
ASTROS.  One  may  also  try  other  elements,  such  as  triangular  or  axisymmetric,  or  other 
nodal  arrangements. 

7.  Results  that  don’t  make  good  physical  sense  should  be  reevaluated.  The  nozzle  stress  results 
are  believable  because  ASTROS  and  MSC/NASTRAN  were  shown  to  return  accurate  values 
and  because  the  odd  results  at  the  throat  do  make  good  physical  sense  under  careful  analysis. 

Thesis  Contribution  And  Suggestions  For  Future  Work 

This  Thesis  has  provided  a  comprehensive  comparison  of  MSC/NASTRAN  and  ASTROS  as 
applied  to  the  analysis  and  optimization  of  shells.  Problems  with  the  porting  of  ASTROS  to  the 
Sun  Workstations  were  revealed.  Finite  element  results  have  been  compared  to  exact  solutions  and 
a  new  nozzle  analysis  has  probed  the  strengths  and  weaknesses  of  each  program. 

The  nozzle  analysis  included  development  of  a  nozzle  model  generation  program.  This  pro¬ 
gram  has  the  flexibility  to  encompass  different  geometries  and  material  properties.  The  analysis 
identified  the  dominant  effects  in  key  regions  of  the  nozzle  and  issues  of  concern  to  the  analyst. 

Work  was  also  done  with  an  axisymmetric  element.  The  potential  for  implementing  this 
element  in  ASTROS  was  demonstrated,  and  a  modularized  FORTRAN  program  developed  for  its 
efficient  generation.  This  flexible  program  can  be  used  alone  to  generate  shell  element  models,  or 
it  could  provide  the  core  for  implementing  axisymmetric  analysis  and  optm  'tion  in  ASTROS. 

Future  work  should  refine  the  nozzle  analysis  to  include  more  accurate  pressure  loads,  com¬ 
posite  materials,  and  a  regenerative  cooling  scheme.  Additional  testing  of  the  axisymmetric  element 
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should  be  done  and  its  incorporation  into  ASTROS  pursued.  Other  areas  of  investigation  could 


include  optimization  with  the  axisymmetric  element  and  its  attachment  to  other  element  types. 


Appendix  a.  caurrmt  iwnox 


Hardware 

Finite  element  analysis  and  optimization  were  performed  on  the  Sun  •!  Workstation  using 
M SC / N ASTRA N  Version  OCA  and  ASTHOS  Versions  -1  and  5.  Models  were  generated  on  tin;  Sun 
•1,  with  tin'  exception  of  the  nozzle  model,  which  required  1MSL  subroutines  available  on  the  VAX. 
Some  post,  processing  used  SDRC/1DFAS  on  tin'  Sun  4. 

Nozzle  Model  Generation  Program 

The  nozzle  generation  program  entailed  significant  effort  and  is  attached  for  reference  pur¬ 
poses.  'I'he  program  inputs  geometry,  material,  thermal,  and  flow  information  and  produces  a  finite 
element  model  of  a  nozzle  using  four  noded  quadrilateral  elements.  'I'he  user  may  opt  to  exploit 
symmetry,  and  determines  whether  an  ASTROS  or  MSC/NASTRAN  data  file  is  generated.  An 
input  file  for  the  axisymmetric  element  program  may  also  he  created,  allowing  the  generation  of  an 
axisymmetric  model  of  the  nozzle. 


PROGRAM  NQGEN 

DIMENSION  ZI (20) ,  RI(20),  SC(20) ,  Y(4),  BREAK(20) ,  CSC0EF(20,4) , 
&ZIK20),  TH(20)  ,  BREAKT(20)  ,  CSC0ET(20,4) 

EXTERNAL  F 
COMMON  GM,  ARATIO 

CHARACTER  INFILE*8,  0UTFILE*8,  TEMP*78,  HEADER*8, 

6EXTN*8 ,  TMP2*8,  LABEL* 1 ,  RTYPE*8  ,  PTYPE*8 

1  F0RMAT(A78) 

2  FORMAT(A) 

C 

C  GET  STARTUP  INFO  FROM  USER 
C 

FP.TNT*.  ’INPUT  FILE  NAME7  ’ 

READ(*  ,2)  INFILE 

PRINT*,  ’OUTPUT  FILE  NAME7’ 

READ(*  ,2)  OUTFILE 
PRINT*,  ’HEADER  FILE  NAME7’ 

READ ( *  ,  2)  HEADER 

PRINT*,  ’1  -  IF  CREATING  A  NASTRAN  DATA  FILE’ 

PRINT*,  ’2  -  IF  CREATING  AN  ASTROS  DATA  FILE’ 

PRINT*,  ’3  -  IF  CREATING  AN  AXYSYMMETRIC  DATA  FILE’ 

READ*,  MODEL 

IF (MODEL . NE . 1 . AND. MODEL  NE .2 . AND. MODEL . NE . 3)  THEN 
PRTNT* ,  ’MUST  HE  ! ,  2  OR  3  "  ’ 

STOP 
END  IF 

IF  (MODEL.  Ft) -3)  GO  TO  !, 

PRINT*,  ’1  -  IF  DESIRE  FULL  MODEL’ 

PR  I  NT  *  ,  ’/  -  IF  DESIRE  QUARTER  MODEL’ 

PRINT*.  IF  DESIRE  STRIP  MODEL’ 


READ* ,  IQUART 

IF  ( I  QUART .  NF. .  1 .  AND .  I  QUART  .  NE .  2  .  AND .  I  QUART .  NE .  3  )  THKN 
PRINT* , 'MUST  BE  1 ,  2  OR  3  !  1  ’ 

STOP 

ENDIF 

r>  PRINT*,  ’MOW  MANY  EI.EMENTS/CIRCUMFERENCE  ??’ 

READ*,  M 

PRINT*,  ’ 1  -  FOR  NO  THERMAL,  2  -  FOR  YES’ 

READ*,  MODTMP 
C 

C  INITIALIZE  FILES 
C 

OPEN (UNIT=6 ,  FILE=INFILE,  STATUS=’OLD’ ) 

OPEN (UNIT=7 ,  FILE=OUTFILE,  STATUS= 'NEW  ’ ) 

OPEN (UNIT=8 ,  FILE=HEADER,  STATUS= 'OLD  ’ ) 

IF(MODEL .EQ .3)  THEN 

OPEN (UNIT=9 ,  FILE=  ’ AXYDATA ’ ,  STATUS3 ’ NEW ’ ) 

ENDIF 

C 

C  INPUT  MODEL  INFORMATION  FROM  INFILE 
C 

READC6 ,  *)  NI  ! NUMBER  OF  INTERPOLATING  POINTS 

READ(6 ,  *)  (ZI(I),  RI(I),  1=1, NI)  (INTERPOLATING  POINTS 

READ(6 ,  *)  NTHRT  ! NTHRT  -  DS  NO  FOR  THROAT 

READ(6 , *)  NIT  (NUMBER  OF  INTERPOLATING  POINTS  FOR  TEMPERATURE 
READ(6,*)  (ZIT(I)  ,  TH(I)  ,  1  =  1, NIT)  ! INTERPOLATING  POINTS 
READC6  ,  *)  TECOEF,  TREE 

READC6 , *)  GM.PO  (READ  IN  GAMMA  AND  STAGNATION  PRESSURE 
READC6  ,  *)  EMOD ,POISS , RHO  (MODULUS,  NU,  AND  MASS  DENSITY 
READC6 , *)  BST,BSC,BSS 

READ(6 , *)  XINIT.XLB.XUB  (DESIGN  VARIABLE  BOUNDS 
READ (6 , *)  XLALL1 .XUALL1  (LOWER  k  UPPER  ALLOWABLES 
READC6 , *)  XLALL2 ,XUALL2 
READC6  ,  *)  XLALL3 , XUALL3 
READ(6 , *)  XLALL4.XUALL4 
C 

C  SET  UP  FOR  QUARTER  MODEL 
C 

IF (IQUART.EQ . 1)  MM=M 
IF( IQUART. EQ . 2)  MM=(M/4)+l 
IF (IQUART.EQ. 3)  MM=2 
C 

C  USE  IMSL  TO  ESTABLISH  SPLINE  FIT 
C 

CALL  CSAKM(NI,ZI,RI, BREAK, CSCOEF)  (SPLINE  FIT  THE  INTERP  PTS 
CALL  CSAKM(NIT,ZIT ,TH .BREAKT .CSCOET)  (FIT  TEMPERATURE 
C 


C  FORMATS 
C 

70  FORMAT ( ’CONROD  ’,4I8,F8.3) 

72  FORMATC 'FORCE  \3I8,4F8.3) 

80  FORMAT! 'DESVARP  ’,2I8,3E8.2) 

81  FORMATC 'PLIST  ’,I8,A8,I8) 

82  FORMATC 'DCONSTR  ’,I8,A8) 

83  FORMATC 'DCONDSP  ’  , 218 , A8  ,  F8 . 3 , A8.2I8 ,F8 . 3) 

90  FORMATC 'EIGR  ’.IS  ,  A8 ,24X ,  18) 

91  FORMATC 'TEMP  > .2I8.F8.1) 

100  FORMATC ’BEGIN  BULK’) 

101  FORMATC 'GRID*  ’ ,21 16 .2E16.9, A3, II ,/ , A3.I1 ,4X ,E16.9 ,2116) 

102  FORMATC 'GRID*  '  , 21 16 ,2E16. 9, A3, II ,/, A3, II ,4X ,E16. 9, 116) 

103  FORMAT ( ’ SPC  ’.3I8.F8.3) 

104  FORMATC 'CQUAD4  ',618) 

100  FORMATC 'PLOAD  ’  , 13 , F8 . 3 , 418) 

106  FORMATC  'MAT1  ’  ,  18 ,  F.8 . 3 , 8X  ,4F.8 . 3,  SX ,  A6,  / ,  AS ,  3E8 . 3) 

107  FORMATC 'PSHELL  ’ , 218 , F8 . 3 , 18 ,8X , 18) 

108  FORMATC 'ENDDATA  ’) 


109 

FORMATC 

DF.SVAR  ’ 

13 

A1  ,11  ,6X, 

3E8.2) 

1  10 

FORMATC 

DVPRFI.l  ’ 

13 

’PSHELL 

’ ,218, 

32X , 

A3 

11 ,/,A3 

T 1 

1  1 1 

FORMATC 

DRESP 1  ’ 

1 3 

A  1 , 1 2 , SX , 

2A8.8X 

,18, 

8X 

16) 

1  12 

FORMATC 

DCONSTR  ’ 

T  H 

'  ALL 

’ ,2F8. 

0) 

1  13 

FORMATC 

DRESP 1  ’ 

I  o 

2A8) 

1  14 

FORMAT ( 

DESOB  I  ’ 

13 

AR.HX.AS) 

1  IS 

FORMAT ( 

DOPTPRM  ’ 

:*  i 

, 2F8 . 3) 

1  16 

FORMATC 

C0RD2C  ’ 

/  i  < 

,  >  i  . .  ,  ’( 

20’  ,/, 

’  +  2( 

'  ,  3F8 

3) 

1  17 

FORMAT ( 

nUPTPRM  ’ 

i.  r 

’.) 

I  1  8 

FORMATC 

PARAM , PUS 

i  ,  - 

(’) 

1  19 

’0  '".MAT ( 

A!  .  I  2 , SX , 

3 EH. 2) 

i  7'; 

(■  '  .RMATC 

'! SHELL 

’.21 K , 

32X 

A  ' 

I.’./.  A3 

I  2 

,F8 . 3) 


121  FORMAT ( ’ DRESP 1  ’,18 , A 1 , 13 ,4X ,2A8 ,8X , 18 ,8X , IS) 

122  FORMAT  ( ’GRID*  ’  , 21 16 , 2H1G  .  9 ,  A3 , 12 ,  /  ,  A3 , 12. , 3X  ,K16 . 9 , 1  16) 

C 

129  FORMAT ( ’DKSVAR  ’ , 13 , A 1  , 13 , 4X , 3K8 . 2) 

130  FORMAT ( ’ DVPREL1  ’  ,  18 , ’PSHEI.L  ’ , 218 , 32X ,  A3 , 13 , / , A3 , 13 , 2X , 18 , FS . 3) 

131  FORMAT ( ’DRESP  1  ’  ,  18  ,  A 1  , 1'l , 3X , 2A8 ,8X  ,  18 ,8X  ,  IS) 

132  FORMAT( ’GRID*  ’ ,2I1G ,2E16. 9, A3, 13 ,/ , A3 , 13 ,2X ,E16. 9 ,1 16) 

C 

139  FORMATC ’DESVAR  ’ , 18, A1 , 14 , 3X.3E8 .2) 

140  FORMAT ( ’ DVPREL1  ’,18,’PSHELL  ’ ,218 ,32X , A3, 14,/ , A3 , 14, IX , 18 ,F8 .3) 

141  FORM AT (’ DRESP 1  ’ , 18 , A1 ,I5,2X,2A8 ,8X ,18 ,8X, 18) 

142  FORMATC ’GRID*  ' , 2116, 2E16. 9, A3, 14 ,/ , A3 , 14 , IX ,E16. 9 ,116) 

C 

1S2  FORMAT ( ’GRID*  ’ ,2116, 2E16. 9, A3.IS ,/, A3.IS.E16.9 ,116) 

C 

C  READ  HEADER  AND  WRITE  TO  DATA  FILE 
C 

IFCMODEL .EQ . 3)  THEN 
NF=9 
ELSE 
NF=7 
ENDIF 

READ(8  ,  *)  NLINES  1INPUT  HEADER  INFORMATION 
DO  155  1=1, NLINES 
READ(8 , 1)  TEMP 

WRITE (NF , 1 )  TEMP  IWRITE  HEADER  TO  FILE 
155  CONTINUE 
C 

C  FIT  MESH  TO  NOZZLE  GEOMETRY 

C 

RM=FLOAT (M) 

PI=3. 141592653589793 
Z=0.0  ! START  AT  END 

NZ=0  INUMBER  OF  ELEMENTS  IN  Z  DIRECTION 
STEP=2. *SIN(PI/RM) 

160  NZ=NZ+1 

R=CSVAL(Z , NI- 1 .BREAK , CSCOEF) 

D=R*STEP  IDEAPTH  TO  STEP 
Z=Z+D 

IF(Z.LT.ZICNI))  GO  TO  160 
IF(Z-ZI(NI) ,LE. D/2.0)  THEN 
FACTOR=ZI(NI)/Z 
ELSE 

FACTOR=ZI(NI)/(Z-D) 

NZ=NZ-1 

ENDIF 

C 

C  START  BULK  DATA  DECK 
C 

WRITECNF, 100) 

IF(MODEL . NE . 3)  THEN 

WRITE (7, 116)  1,0,0.,0.,0.,0.,0.,1.,1.,0.,1. 

ENDIF 

C 

C  MATERIAL 
C 

WRITE(NF  ,  106)  1 , EMOD .POISS  ,RHO .TECOEF.TREF, ’ +MT1  ’, 
ft’+MTl  ’ ,BST,BSC,BSS 

C 

C  FIRST  GRID  RING  WITH  SPC’S 
C 

IF (MODEL. EQ. 3)  THEN 

WRITE (9, 102)  1 , 0 , RI ( 1 ) ,0.  ,  ’+GR’ ,1 , ’*GR’ ,1 ,0.  ,0 
WRITEO ,  103)  1,1,123456,0. 

ELSE 

DO  165  1=1 , MM 

PX=RI(1)*C0S(2. *P I* FLO AT ( I —  1 ) / RM ) 
PY=RI(1)*SIN(2.*PI*FI.0AT(I-1)/RM) 

IF (I .LT. 10)  THEN 

WRITE  (7,  102)  1 , 0 , 1’X  ,  PY  ,  ’+GR  ’ ,  I  ,  ’  *GR  ’,1,0.  ,  1 
ELSEIFCI . LT. 100)  THEN 

WRITE (7, 122)  1 ,0 ,PX  ,PY  ,  ’ +  GR ’ , I , ’*GR’ ,1,0.  ,  1 
EI.SF.TE ( I  .  I.T .  1000)  THEN 

WRITE (7 , 1 32)  I.O.PX.PY,  ’ +GR ’ , I , ’ *GR ’ . I  ,0.  ,  1 
ELSEIFCI  .  I.T.  10000)  THEN 

WRITE  (7 .142)  I  ,0,1’X  ,PY,  ’  *  GR  ’ ,  I  ,  ’  *GR  ’  ,  I  ....  1 
ELSEIFCI . GE . 10000)  THEN 

PRINT*.  ’LOGIC  ERROR  LOOP  165’ 


STOP 

F.NOIF 

WRITE  (7 ,103)  1  , 1, 23456,0.  !  WRITE  SPC’S 
I F (MODTMP  .  HQ . 2)  THEN 
WRITE (7 ,91)  5 , 1 ,TH (  1 ) 

ENDIF 

165  CONTINUE 
ENDIF 
C 

C  WRITE  AXYELEM  FIRST  DATA 
C 

IF (MODEL . NE. 3)  GO  TO  167 

WRITE (7 , *)  DBLE(EMOD) .DBLE(POISS) 

WRITE (7,*)  NZ 

WRITE(7  ,  *)  DBLE(ZKl))  ,DBLE(RI(1)) 

167  CONTINUE 
C 

C  INITIALIZE  VALUES  FOR  LOOPING  TO  END  OF  NOZZLE 
C 

ATHRT=RI (NTHRT) »RI (NTHRT)  ! THROAT  AREA  WITHOUT  PI 
Z2=0 . 

NZA=0  'THIS  IS  THE  NUMBER  OF  TIMES  STEPPED  IN  Z  DIRECTION 
D=RI (1) *STEP 
C 

C  LOOP  TO  END  OF  NOZZLE 
C 

170  NZA=NZA+1 
Z1=Z2 
Z2=Z2+D 
Z1C=Z1*FACT0R 
Z2C=Z2*FACT0R 

R=CSVAL(Z2,NI-1 .BREAK .CSCOEF) 

THETMP=CSVAL (Z2C , NIT- 1 , BREAKT , CSCOET) 

D=R*STEP  ! SET  D  FOR  NEXT  LOOP 
R2C=CSVAL(Z2C,NI-1, BREAK .CSCOEF) 

IF (R2C.I.E.RI (NTHRT))  THEN 
R2C=RI (NTHRT) 

ENDIF 

IF(M0DEL.EQ.3)  THEN 
II=NZA*M+1 
IF(II.LT.IO)  THEN 

WRITE(9 , 102)  II.0.R2C.0. ,,+GR,.II, ’ *GR  \  II ,Z2C , 0 
ELSEIF(II.LT.IOO)  THEN 

WRITE(9 , 122)  II.0.R2C.0. . '+GR’,II , ’ ♦GR'.II ,Z2C .0 
ELSEIFdl.LT.  1000)  THEN 

WRITE(9 , 132)  II.0.R2C.0. , ’ +GR’ ,11, ’*GR’ .II.Z2C.0 
ELSEIFdl.LT.  10000)  THEN 

WRITE(9 , 142)  II ,0 ,R2C ,0. , ’+GR’ , II , ' *GR’ , II ,Z2C,0 
ELSE 

PRINT* , 1  LOOP  TOO  BIG  AFTER  170’ 

STOP 

ENDIF 

WRITE(9 , 103)  1,11,246,0. 

TI=NZA*M-M+1 

WRITE(9 ,70)  II , II , II+M, 1,1.  ! DUMMY  CONROD 
ELSE 

DO  175  1=1, MM 

C  PX=R2C*C0S(2 . *PI*FLOAT (I-1)/RM) 

C  PY=R2C*SIN(2 . *P I* FLO AT (I-1)/RM) 

THETA=FLOAT(I- 1) *360 .  /RM 

IC0RD=1 

II=NZA*M+I 

IFdl.I.T.  10)  THEN  IWRITE  GRID  RING 

WRITE(7 , 102)  1 1 , ICORD , R2C .THETA , ’+GR’ , II , ’ *GR ’  ,  II , Z2C , I  CORD 
ELSEIFdl.LT.  100)  THEN 

WRITE (7, 122)  II .ICORD, R2C, THETA, ’*GR ’ ,11. ’ *GR ’ , II ,Z2C , ICORD 
ELSEIFdl.LT.  1000)  THEN 

WRITE (7, 132)  IT .ICORD, R2C, THETA, ’ +GR ’ . II,  *GR’ .II.Z2C , ICORD 
ELSEIFdl  . I.T.  10000)  THEN 

WRITE  (7, 142)  1 1,  ICORD,  R2C,  THETA,  ’-*GR’  .11  ,  '  *GR  ’  ,  1 1 ,  Z2C  ,  ICORD 
ELSEIFdl  .I.T.  100000)  THEN 

WRITE  (7  ,  152)  II ,  ICORD,  R2C.  TIIF.TA,  ’ +GR’  ,  II  .  ’  *GR  ’  ,11  ,Z2G  .ICORD 
ELSEIFdl  .GE.  100000)  THEN 

PRINT*.  ’LOGIC  ERROR  LOOP  175’ 

STOP 
END  IF 

I  F  ( MnDTMi’ .  HU  ■  2)  THEN 

WR I TF (7,0l)  5,11 , THETMP 
ENDIF 


IF( I  QUART . EQ. 2)  THEN  'WRITE  SEC'S  IF  QUARTERED  MODEL 
1FCI.EQ.1)  THEN 

WRITE(7 , 103)  1,11,246,0.0 
KLSEIF(I.EQ.MM)  THEN 

WRITE (7 ,103)  1,11,246,0.0 
ENDIF 

ELSE  IF ( I QUART . EQ . 3)  THEN  ! WRITE  SEC’S  FOR  STRIP  MODEL 
WRITEC7 , 103)  1,11,246,0. 

ENDIF 

17:  CONTINUE 

DO  177  1  =  1, MM- 1  ! WRITE  QUAD  ELEMENTS 

II=NZA*M-M+I 

WRITEC7, 104)  II, NZA, 11,11+1 , II+M+1 , II+M 

177  CONTINUE 

IF(IQUART.EQ.l)  THEN  ! IF  NOT  QUARTERED,  COMPLETE  RING 
WRITE(7 , 104)  NZA«M, NZA , NZA*M , NZA+M-M+ 1 , NZA»M+ 1 , NZA*M+M 
ENDIF 

WRITE(7 , 107)  NZA,  1,1.,  1,1  ( WRITE  PSHELL  FOR  RING 

ENDIF 

178  ZP=(Z2C+ZlC)/2.0  (Z  LOCATION  FOR  PRESSURE 
RP=CSVAL(ZP , NI-1 , BREAK, CSCOEF)  (RADIUS  AT  Z  LOCATION 
IF(RP.LE.RKNTHRT) )  THEN  ! LONT  LET  RP<GIVEN  R  THROAT 

RP=RI (NTHRT) 

ENDIF 

ARATIO=RP»RP/ATHRT  (AREA  RAIIO 
PP0T=(2 . O/CGM+1 .0) )*» (GM/ (GM-1 . 0) ) 

IF(ZP.EQ.ZI( NTHRT ) )  THEN  'ITERATE  TO  PRESSURE 
P=PO*PPOT 

ELSElFCZP.LT. ZI (NTHRT))  THEN 
A=PPOT 
B=1 .0 

135  IF(F((B+A)/2.0) ,LE. 0.0)  THEN 

B=(B+A)/2.0 
ELSE 

A=(B+A)/2.0 

ENDIF 

IF(CB-A) ,GE. 0.00001)  GO  TO  135 
P=PO* (B+A) /2 . 0 
ELSE 
A=0.0 
B=P?OT 

137  IF(F((A+B)/2.0) .LE.O.O)  THEN 

B=(A+B)/2.0 
ELSE 

A=(A+B)/2.0 

ENDIF 

IF(CB-A) .GE. 0.00001)  GO  TO  137 
P=PO*(A+B)/2.0 
ENDIF 

IFCM0DEL.EQ.3)  THEN 

WRITE (7 , *)  DBLECZ2C) ,DBLE(R2C) ,DBLE(P) 

WRITEC9 ,72)  1,NZA*M+1,0,1.  ,1.  ,1.  ,  1 .  (DUMMY  LOAD 
GO  TO  210 
ENDIF 

DO  190  1=1 ,MM-1  (WRITE  PLOAD  ON  ELEMENTS 
II=NZA+M-M+I 

WRITEC7 , 105)  1 ,P, 11,11+1,  II+M+1, II+M 
190  CONTINUE 

IF (IQUART.EQ . 1)  THEN  (IF  NOT  QUARTERED  MODEL,  COMPLETE  PLOAD 
WRITEC7, 105)  1 ,P,NZA*M,NZA*M-M+1 ,NZA*M+1 ,NZA*M+M 
ENDIF 

IFCMODEL.EQ.  1)  THEN  (WRITE  DESVAR,  DVPRF.Ll  FOR  NASTRAN 
IFCNZA.LT. 10)  THEN 

WRITEC7 , 109)  NZA, >T’ , NZA ,X1NIT, XLB ,XUB 
WRITE (7 ,110)  NZA  ,  NZA  ,4 ,  ’+DP  ’  ,  NZA ,  ’  +  DP  ’  ,  NZA  ,N7.A  ,  1  . 
ELSEIFCNZA . LT. 1 00)  THEM 

WRITE (7, 1 19)  NZA, ’T’ , NZA , XI HIT , XLB , XUB 
WRITE (7 , 120)  NZA, NZA, 4, ’ +  DP ' ,NZA, ’ +DP ’ , NZA , NZA ,  1  . 
ELSEIECNZA.LT. 1000)  THEN 

WRITE (7 ,129)  NZA, >T’ , NZA , X  IN  IT , XLB , XUB 
WRITE  (7 , 120)  NZA  ,  NZA  ,4  ,  ’  +DP  ’  ,NZA  ,  ’  +1)1”  ,  NZA  ,  NZA,  1  . 
EI.SKIFCNZA.lt.  10+00)  THEN 

WRITE! 7 , 1  30  j  NZA  ,  ’ T  ’  ,  NZA  ,  X  MIT,  XLB  ,  XUB 
WRITE (7  ,  14  +  )  NZA, NZA, 4, ’  +  DI”  .NZA. '  +  DP  ’  , NZA .NZA, 1  . 
ELSEIFCNZA  .GE.  100';")  THEN 

PRINT*,  ’  LOGIC  ERROR  BETWEEN  1""  AND  2" 

STOP 
END  I F 


O  O  U  OUUCN 


N0S=2  !  NUMBER  OK  STRESS  COMPONENTS  FOR  CONSTRAINTS 

DO  200  1=1, NOS  ! WRITE  DRESP1  FOR  NASTRAN 
IF(I.EQ.l)  THEN 
I.ABEL=  ’  S  ’ 

RTYPE=’  STRESS  ’ 

PTYPE= ’PSHELL  ’ 

IATTA=9 

IATT1=NZA 

ELSEIFCI.EQ.2)  THEN 

labf:l=>s’ 

RTYPE=’ STRESS  ' 

PTYPE=’ PSHELL  ’ 

IATTA=17 

IATT1=NZA 

ENDIF 

IF( ( 10*NZA+I) . LT. 100)  THEN 

WRITEC7 ,111)  NOS*NZA-NOS+I ,  LABEL ,  10*NZA+I , RTYPE.PTYPE , 
ft  IATTA  IATT1 

ELSEIFC ( 10*NZA+1) . LT . 1000)  THEN 
WRITE (7, 121)  NOS*NZA-NOS+I, LABEL, 10*NZA+I, RTYPE.PTYPE, 
ft  IATTA  IATT1 

ELSEIFC (10*NZA+1) .LT. 10000)  THEN 
WRITEC7, 131)  NOS*NZA-NOS+I, LABEL, 10*NZA+I, RTYPE.PTYPE, 
ft  IATTA  IATT1 

ELSEIFC (10*NZA+1) .LT. 100000)  THEN 
WRITEC7, 141)  N0S*NZA-N0S+I, LABEL, 10*NZA+I, RTYPE.PTYPE, 
ft  IATTA, IATT1 
ENDIF 

200  CONTINUE 

WRITE C7 ,112)  N0S*NZA-1, XLALL1 ,XUALL1  (WRITE  DCONSTR  FOR  NASTRAN 
WRITE C7, 112)  NOS+NZA ,XLALL2, XUALL2 
ELSEIF CMODEL . EQ . 2)  THEN  (WRITE  DESVARP ,  PLIST  FOR  ASTROS 
WRITEC7 , 80)  NZA ,NZA ,XLB ,XUB ,XINIT 
WRITEC7, 81)  NZA, ’PSHELL  ’ ,NZA 
ENDIF 

210  IFCNZA.LT.NZ)  GO  TO  170  (LOOP  BACK  IF  NOT  TO  END  OF  MODEL 

NOW  THAT  AT  END  OF  NOZZLE,  FINISH  OFF  MODEL 

IFCMODEL.EQ. 1)  THEN  (FINAL  NASTRAN  DRESP1,  DESOBJ,  DOPTPRM,  PARAM 
PTYPE= ’ 

WRITEC7 ,  111)  88888, ’D’, 88, ’DISP  ’ ,PTYPE, 1 ,NZ*M+1 
WRITEC7 ,111)  99999,  ’D\ 99,  ’DISP  ’  ,PTYPE,  3  ,NZ*M+1 
WRITEC7 , 112)  88888, XLALL3 , XUALL3 
WRITEC7 .112)  99999, XLALL4 , XUALL4 
WRITEC7.113)  N0S*NZ+1,’W  ’.’WEIGHT  ’ 

WRITEC7 , 114)  N0S*NZ+1,’W  ’ , ’MIN 

WRITEC7 ,115)  2,3,20, .5, .01 
WRITEC7 , 118) 

ELSEIF (MODEL. EQ. 2)  THEN  (WRITE  FINAL  DCONSTR  FOR  ASTROS 
WRITEC7 ,82)  1,’VMISES  ’ 

WRITEC7, 83)  1,1, ’LOWER  ’,XLALL3,’D1  ’  ,NZ*M+1 , 1 , 1 . 

WRITEC7, 83)  1,2, ’UPPER  ’,XUALL3,’D2  ’  ,NZ*M+1  , 1 , 1  . 

WRITEC7, 83)  1,3, ’LOWER  ’,XLALL4,’D3  ’ ,NZ«M+1 ,3, 1. 

WRITEC7, 83)  1,4  ’UPPER  ’,XUALL4,’D4  ’ ,NZ*M+1 , 3 , 1 . 

ELSEIF (MODEL . EQ . 3)  THEN 
GO  TO  250 
ENDIF 

WRITF.(7,90)  75,  'GIV  ’,5  (EIGR  INSTRUCTIONS 
WRITE (7 , 108)  (ENDDATA  STATEMENT 

NOW  CLOSE  OFF  FILES  AND  EXIT  PROGRAM 

50  CL0SEC6) 

CL0SEC7) 

Cl.OSE  (8) 

IFCMODEL.EQ. 3)  THEN 
CL0SEC9) 

ENDIF 

STOP 

END 

C  IDEAL  CONVERGENT/DrVF.RGF.NT  NOZZI.F.  EQH 

v, 

REAL  FUNCTION  FIX) 

COMMON  CM,  A  HAT  1(1 
A  CM-t  .0 
!'  CM*  !  .<> 

<:  H/A 


D=A/GM 

E=2.0/GM 

F=((2./B)*»C*A/(X**E*(1 .-X**D)))-ARATIU*ARATIO 

RETURN 

END 


Axisymmrtnc  Finite  Fitment  Generation  Program 

The  axisymmotric  element  program  also  involved  significant  efrort  and  is  attached  for  ref¬ 
erence.  The  FORTRAN  program  inputs  geometry  and  material  properties  and  produces  global 
stiffness,  load  and  stress  recovery  matrices  in  direct  matrix  form.  It,  can  he  easily  modified  since 
it  is  modularized  and  well  documented  internally.  A  Gaussian  quadrature  integration  scheme  is 
used,  and  the  order  of  integration  and  weighting  parameters  can  he  easily  changed  in  the  preface. 
Double  precision  is  used  throughout  to  retain  accuracy. 


PROGRAM  AXYSYM 

DOUBLE  PRECISION  E ,T( 100) . DNU, X ( 100) ,R( 100) ,PR( 100) ,DKGELM(6 ,6) , 
ftPWT(2 , 10) , TEMP (21) ,TEMP2(6) ,TEMP3(14,9) ,SREC(7,6) , 
ftTEMP4(9,3) 

CHARACTER  IHFILE*8 ,0UTFILE*8 
C 

C  EXPLANATION  OF  VARIABLES 
C  E  ==  MGUUI.US 

C  T  ==  VECTOR  CONTAINING  ELEMENT  THICKNESSES 

C  DNU  ==  POISSONS  RATIO 

C  Xi,  Ri  ==  VECTORS  CONTAINING  COORDINATES  OF  END  POINTS 
C  PR  ==  VECTOR  CONTAINING  PRESSURE  ON  ELEMENTS  NORMAL  &  CENTERED 
C  DKGELM  ==  ELEMENT  STIFFNESS  MATRIX,  GLOBAL 
C  NOD  ==  NODE  CONNECTIVITY  TABLE 
C  NG  ==  NUMBER  OF  GRID  POINTS 
C  NE  ==  NUMBER  OF  ELEMENTS 

C  N0Q  ==  ORDER  OF  GAUSSIAN  QUADRATURE 

C  PWT  ==  ARRAY  CONTAINING  PSI  AND  WEIGHTS  FOR  QUADRATURE  (I) 

C  TEMPi  ==  TEMPORARY  STORAGE 

C  PELM  ==  MATRIX  W/C0LUMNS  BEING  ELEMENT  LOAD  VECTORS,  GLOBAL 

C  SREC  ==  STRESS  RECOVERY  MATRIX 

C 

C  ESTABLISH  GAUSSIAN  QUADRATURE  DATA 
C 

N0Q=7 

PWT ( 1 , l)=-0. 949107912342758524526189684048 
PWT(2 , 1)=0 . 1 29434966 1G886969327061 1432679 
PWT(1 ,2)=-0. 741531 185599394439863864773281 
PWT(2, 2) =0.27970539148927606790 146777 1424 
PWT ( 1 , 3) =-0 . 405845 1 5 1 377397 1 669066064 1 2077 
PWT (2 ,3) =0.381330050505118944900369775489 
PWT ( 1 ,4) =0 . OD+OO 

PWT ( 2 , 4 ) =0 . 4 1 7959 1 83673469387755 1 020408 1 6 
PWT (1,5) =-PWT ( 1 ,3) 

PWT (2 , 5) =PVT(2 ,3) 

PWT ( 1  , 6) =-PWT ( 1 ,2) 

PWT (2 , 6) =PWT (2,2) 

PWT ( 1 , 7 ) = -PWT ( 1 , t) 

PWT (2,7) -PWT (2,1) 

<:  GET  STARTUP  PILE  INFORMATION  FROM  USER 
<: 

2  FORMAT f  A ) 


4  F0RHAKE16.8) 

PRINT*, ’INPUT  FILE  NAME  ??’ 

READ (*  ,2)  INFILE 

PRINT*. ’OUTPUT  FILE  NAME  ?? > 

READ ( * , 2)  QUTFILE 
C 

C  INITIALIZE  FILES 
C 

OPEN (UNIT=6,  FILE=INFILE,  STATUS= ’OLD ’ ) 
OPEN (UNIT=7 ,  FILE=OUTFILE,  STATUS= ’ NEW ’ ) 
C 

C  INPUT  DATA 
C 

READ(6  ,  *)  E , DNU 
READC6 ,  *)  NE 
READC6 , *)  X ( 1) ,R(1) 

DO  10  1=1 ,NE 

READ(6 , *)  X(I+1),R(I+1),PR(I) 

10  CONTINUE 

DO  20  1=1, NE 
READ(6 , *)  T(I) 

20  CONTINUE 

C 

C  CREATE  AND  OUTPUT  KGG 


OPEN (UNIT=8 ,  FILE=’ SCRATCH 1> ,  STATUS= ’ NEW ’ ) 

DO  200  1=1 ,NE+1 
IF(I.EQ. 1)  THEN 

CALL  KELM(E,T(I) ,DNU,X(I) ,R(I) ,X(I+1) ,R(I+1) ,NOQ,PWT .DKGELM) 

DO  150  J=1 ,3 
WRITE(8,3)  2 
WRITE(8 , 3)  2* J-l 
WRITE (8 ,3)  2 
WRITE (8 ,3)  1 
DO  148  K=1 ,6 
WRITE(8 ,3)  1 
WRITEC8 ,4)  DKGELM(K.J) 

WRITEC8 ,3)  1 
WRITE (8 ,4)  O.OD+OO 

148  CONTINUE 
IF(NE.EQ.l)  GO  TO  149 
DO  149  K=13 , NE*6+6 

WRITE(8 ,3)  1 
WRITE(8 ,4)  O.OD+OO 

149  CONTINUE 
WRITEC8 , 3)  2 
WRITE(8 ,3)  2* J 
WRITEC8 ,3)  2 
WRITE(8 ,3)  1 

DO  150  K=1 , NE*6+6 
WRITEC8 ,3)  1 
WRITE (8 ,4)0. OD+OO 

150  CONTINUE 
ELSEIF(I.EQ.NE+1)  THEN 

CALL  KELM(E,T(I-1), DNU, Xd-1),R(I-1),X(I),R(I),N0Q,PWT, DKGELM) 
DO  160  J=4 , 6 
WRITE(8 ,3)  2 

WRITE (8 ,3)  2* (3*NE) -5+2*( J-l) 

WRITE (8 ,3)  2 
WRITEC8 ,3)  1 
IFCNE.EQ. 1)  GO  TO  154 
DO  154  K=1 ,6*NE-6 
WRITE (8, 3)  1 
WRITE(8 ,4)  O.OD+OO 
154  CONTINUE 

DO  158  K=  1  , 6 
WRITE (8 , 3)  1 
WRITE (8, 4)  DKGELM (K,J) 

WRITE!?., 3)  1 
WRITE! 8 ,4)  O.OD+OO 
158  CONTINUE 

WRITE (8,3)  2 
WR I TE (8,3)  6*NF.+2*  1-6 
WRIT  Ed.  3)  2 
WRITE (8.0  1 
IK!  16!)  K"  I  ,  NE*6  *6 
WRITE! 8,8;  i 
WHITE!  8, 41 
166  CONTINUE 


165 

166 


168 

170 


178 

180 


184 


188 


189 


190 

200 


ELSK 

DO  166  J=1 ,3 
DO  165  K  =  1 , 9 

TEMP4 (K  ,  J)=0 . OD+OO 
CONTINUE 


CONTINUE 

CALL  KELM(E,T(I-1),DNU,X(I-1)  ,R(I- 1)  ,  X(I)  ,K(I)  ,  NOQ  .PUT,  DKGEI.M) 
DO  170  9=1,3 
DO  168  K=1 , 6 

TEMP4(K,  J)=DKGELM(K, J+3) 

CONTINUE 

CONTINUE 

CALL  KELM(E,T(I) ,DNU,X(I) ,R(I) ,X(I+l) ,R(I+1) , NOQ ,PWT .DKGELM) 

DO  180  J=1 ,3 
DO  178  K=1 , 6 

TEMP4 ( K+3 , J) =TEKP4 (K+3 , J)+DKGELM(K  ,  J ) 


CONTINUE 
CONTINUE 
DO  190  J=1 ,3 
WRITE(8 ,3)  2 

WRITEC8 ,3)  2*(3*I)-5+2*(J-l) 

WRITEC8 ,3)  2 
WRITE(8,3)  1 
IFCI.EQ.2)  GO  TO  184 
DO  184  K=l, 6*1-12 
WRITE (8 ,3)  1 
WRITEC8 ,4)  O.OD+OO 
CONTINUE 
DO  188  K=1 , 9 
WRITE(8 ,3)  1 
WRITEC8  ,4)  TEHP4(K , J) 

WRITEC8 ,3)  1 
WRITE(8 ,4)  O.OD+OO 
CONTINUE 

IF(I.Eq.NE)  GO  TO  189 
DO  189  K=6*I+8,NE*6+6 
WRITE (8 ,3)  1 
WRITEI8 ,4)  O.OD+OO 
CONTINUE 
WRITEC8 , 3)  2 
WRITEC8 , 3)  6*1+2* J-6 
WRITEC8 ,3)  2 
WRITEC8 ,3)  1 
DO  190  K=1 , NE*6+6 
WRITE(8 ,3)  1 
WRITE(8,4)  O.OD+OO 
CONTINUE 
ENDIF 
CONTINUE 
WRITEI8 , 3)  3 

CALL  ASTOUTC 'HYKGG  ’,’RDP  ’,’REC  ’ ,NE*6+6 , NE*6+6 , 
ft ’SCRATCH1 ’  ,7) 

CL0SE(8) 


C 

C  CREATE  NODAL  LOAD  VECTOR 
C 

OPEN (UNIT=8 ,  FILE= ’SCRATCH2 ’  ,  STATUS* ’ NEW ’ ) 

0  300  1=1 ,NE+1 
IFCI.EQ. 1)  THEN 

CALL  PL'JAD(X(I)  ,R(I)  ,X(I  +  1)  ,R(I  +  1)  ,PR(I)  ,NOq  .PWT.TEMP2) 
WRITEI8 , 3)  2 
WRITEC8  3)  1 
WRITEI8 , 3)  2 
WRITEI8 , 3)  1 
DO  220  J=1 ,3 
WRITE (8, 3)  1 
WRITE (8, 4)  TF.MP2 ( J ) 

WRITEC8 ,3)  1 
WRITE (8 , 4 )  O.OD+OO 
220  CONTINUE 

KL5KI F ( I . EQ . NE+ 1 )  THEN 

CALL  PL0AD(X(I-1)  ,  R(I-l)  ,X(I)  ,R(I)  ,PR(I-1)  .  NOq  ,  PWT  .TF.MP2) 


DO  230  J  =  4 , 6 
WR ITE (8 , 3)  1 
WRITE (8 , 4)  I EMP2 ( .1) 
WRITE (8, 3)  1 
WR I TK (8,41  O.OD+OO 
CONTINUE 


ELSE 


230 


o  o  o 


CALL  PLOAD(Xd-l)  ,R(I-1)  ,X(I)  ,R(I)  ,PRd-l)  ,  NOQ  ,PWT ,TEMP2) 
TEMP ( 1 ) =TEMP2 (4) 

TEMP  (  2  )  =TKMP2  ( 5>) 

TEMP ( 3 ) =TEMP2 ( C) 

CALL  PLOAD(X(I)  ,R(I)  ,  X  ( I  + 1 )  ,R(I  +  1),PR(I)  .N0Q.PWT.TKMP2) 
TEMP ( 1 )=TEMP( 1 ) +TEMP2 ( 1 ) 

TEMP (2) =TEMP ( 2) +TEMP2 (2) 

TEMP(3)=TEMP(3)+TEMP2(3) 

DO  240  J  =  1 ,3 
WRITE(8,3)  1 
WRITE(8 ,4)  TEMP ( J) 

WRITE (8 , 3)  1 
WRITEC8.4)  O.OD+OO 
240  CONTINUE 

ENDIF 

300  CONTINUE 

URITE(8,3)  3 

CALL  ASTOUTC’MYPL  ’  ,’RDP  ’ ,’REC  \NE*C+6,1, 
ft ’SCRATCH2 ’ ,7) 

CL0SE(8) 

CREATE  AND  OUTPUT  STRESS  RECOVERY  MATRIX 

OPEN(UNIT=S,  FILE= 'SCRATCH3' ,  STATUS= 'NEW') 

DO  400  1=1 ,NE+1 
IF(I.EQ.l)  THEN 

CALL  RECOVER(X(I)  ,R(I)  ,X(I+1)  ,R(I+1)  ,DNU,E,T(I)  ,SREC) 

DO  330  J=1 ,3 
WRITEC8 ,3)  2 
WRTTF. (8 , 3)  9*3-1 
WRITEC8 , 3)  2 
WRITEC8.3)  1 
DO  320  K=1 ,7 
WRITE  (8, 3)  1 
WRITE(8,4)  SREC(K.J) 

32C  CONTINUE 

IF(NE.EQ.l)  THEN 
GO  TO  329 
FI  SF 

DO  328  K=8,7*NE 
WRITE (8, 3)  1 
WRITE (8 ,4)  O.OD+OO 
3  ,  CONTINUE 

ENDIF 

3"0  CONTINUE 

WRITE (8, 3)  2 
WRITE(8 ,3)  2* J 
WRITE(8 ,3)  2 
WRITEC8 ,3)  1 
DO  330  K=1,7*NE 
WRITE (8 ,3)  1 
WRITEC8 ,4)  O.OD+OO 
3  )  CONTINUE 

ELSEIFCI .EQ . NE+1 )  THEN 

CALL  RECOVER(XCI-l) ,R(I-1) ,X(I) ,R(I) ,DNU,E,Td-l ) ,SREC) 
DO  340  J=4 , 6 
WRITE(S,3)  2 

WRITE(8 , 3)  l+6«NF.+2*(J-4) 

WRITEC8 , 3)  2 
WRITE(8 , 3)  1 
IF(NE . EQ .  1)  GO  TO  334 
DO  334  K=  1 , 7* (NE- 1) 

WRITE (8 ,3)  1 
WRITE(8 ,4)  O.OD+OO 
3..  I  CONTINUE 

DO  338  K  =  1 ,7 
WRITE (8 ,3)  1 
WRITE (8, 4)  SREC(K.J) 

338  CONTINUE 

WRITE (8 , 3)  2 
WRITE(8 , 3)  G*NE+2*.I-G 
WRITE (8. 3)  2 
WRITE (8. 3)  1 
DO  .140  K '  1  .  7  •  N E 
WRITE 7 8, 3)  I 
WRITE  18, 4)  O.OD  +  OO 
CONTINUE 
ELSE 

CALL  RECOVER!  Xd  -  U  , R(  i -  1 )  .  X ( I )  , kd  )  . i  N- : . E , T d  1  )  .ELEC) 


DO  350  K  =  1 , 7 
DO  3 '18  J  =  l,6 

TF.MP3  (K ,  J )  =SREC (K  ,  J) 

348  CONTINUE 

350  CONTINUE 

CAI.L  RECOVER(X(I)  ,R(I)  ,X(I  +  1)  , R(  1  +  1 )  , DNU,E,T( I )  ,SREC) 

DO  360  J=  1 , 3 
WRITE(8 ,3)  2 

WRITEC8 ,3)  6*(I-l)+2* J-l 
WRITE(8 ,3)  2 
WRITE(8 ,3)  1 
IF(I.EQ.2)  GO  TO  354 
DO  354  K=l, 7*1-14 
WRITE(8 ,3)  1 
WRITE(8,4)  O.OD+OO 
354  CONTINUE 

DO  356  K=l,7 
WRITE (8, 3)  1 
WRITE(8 ,4)  TEMP3 (K , J+3) 

356  CONTINUE 

DO  358  K=1 ,7 
WRITE(8,3)  1 
WRITEC8 ,4)  SREC(K , J) 

358  CONTINUE 
IF(I.EQ.NE)  GO  TO  359 
DO  359  K=7*I+2,7*NE 

WRITEC8 ,3)  1 
WRITE(8 ,4)  O.OD+OO 

359  CONTINUE 
WRITE (8 , 3)  2 
WRITE(8,3)  6*1+2* J-6 
WRITE(8 , 3)  2 
WRITE! 8 , 3)  1 

DO  360  K=1,7*NE 
WRITE! 8, 3)  1 
WRITEC8 ,4)  O.OD+OO 

360  CONTINUE 
ENDIF 

400  CONTINUE 

WRITE(8,3)  3 

CALL  ASTOUTC ’MYSREC  ',’RDP  ’ , ’REC  ’  , 7*NE ,NE*6+6 , 
ft ’ SCRATCH3 ’ ,7) 

CL0SE(8) 

STOP 

END 

C******************** ****************************************************** 

C  THIS  SUBROUTINE  INPUTS  GEOMETRY  AND  MATERIAL  DATA  AND  RETURNS  AN  ELEMENT 
C  STIFFNESS  MATRIX  (DKELM)  IN  GLOBAL  COORDINATES. 

C 

SUBROUTINE  KELM(E,T,DNU,X1 ,R1,X2.R2,N0Q ,PWT,DKGELM) 

DOUBLE  PRECISION  E.T.DNU, XI ,R1 , X2 , R2 ,PWT(2, NOQ)  , 
ftDK(6 ,6) ,DKGELM(6,6) ,B(5,6) ,G(5,5) ,DL 
C 

C  EXPLANATION  OF  VARIABLES,  (I)  INPUT,  (0)  OUTPUT 
C  E  ==  MODULUS  (I) 

C  T  ==  ELEMENT  THICKNESS  (I) 

C  DNU  ==  POISSONS  RATIO  (I) 

C  Xi,  Ri  ==  COORDINATES  OF  END  POINTS  (I) 

C  NOQ  ==  ORDER  OF  GAUSSIAN  QUADRATURE  (I) 

C  PWT  ==  ARRAY  CONTAINING  PSI  AND  WEIGHTS  FOR  QUADRATURE  (I) 

C  DKGELM  ==  ELEMENT  MATRIX,  GLOBAL  (0) 

C  DK  ==  ELEMENT  MATRIX  FOR  INTERNAL  USE,  SQUARE,  LOCAL  COORDINATES 

C  B  ==  B  MATRIX 

C  G  ==  E  MATRIX 

C  DL  ==  LENGTH  OF  ELEMENT 

C 

C  CLEAR  DK  MATRIX 
C 

DO  10  1=1,6 
DO  8  1=1,6 

DK(I  ,.J)  =0.00+00 
8  CONTINUE 

10  CONTINUE 

C 

0  PERFORM  GAUSSIAN  QUADRATURE  TO  FORM  ELEMENT  K  MATRIX  ,  DI./2  IS  1 

C 

DO  100  1=1, NOQ 

CALL  RMAT (  X  1  ,kl  ,  X2  ,  R2  , PWT <  1,1)  ,B,I(L) 

CALL  EM AT ( DNU , E ,T ,G) 


II 


oooooooooooooo 


DO  90  J 1=  1 ,  G 
DO  90  J2=  1 , 5 
DO  90  J3= 1  ,5 
DO  90  J4  =  1  ,  G 

DK(J1 , J4)=DK(J1 , J4)+PWT(2,I)*B(J2, J1)*G(J2, J3)*B(J3, J4) 

90  CONTINUE 

100  CONTINUE 

DO  lbO  1=1,6 
DO  145  J=1 ,G 

DK  ( I ,  J)  =DK (I ,  J) *DL/2 .OD+OO 
145  CONTINUE 

150  CONTINUE 
C 

C  ROTATE  TO  GLOBAL  COORDINATES 
C 

CALL  ROTATE (X 1 , R1 , X2 , R2 , 6 , 6 , DK ,  5.DKGELM) 

C 

C  FINISHED 
C 

RETURN 

END 

C************************«************************************************* 

C  THIS  SUBROUTINE  INPUTS  GEOMETRY,  PRESSURE,  AND  GAUSSIA*'  DIRECTIONS,  AND 
C  OUTPUTS  ELEMENTAL  NODAL  LOAD  VECTOR  IN  GLOBAL  COORDINAii  3. 

C 

SUBROUTINE  PLOADCX 1 , R1 , X2 , R2 ,PR ,NOQ ,PWT ,PGENL) 

DOUBLE  PRECISION  XI  ,R1  ,X2  ,R2  ,PR,PtfT(2,N0Q)  .PGENI/S)  ,PEENL(6)  , 
ftDN (3 ,6)  ,DL,P(3) 

EXPLANATION  OF  VARIABLES,  (I)  INPUT,  (0)  OUTPUT 
Xi,  Ri  ==  COORDINATES  OF  END  POINTS  (I) 

PR  ==  STATIC  PRESSURE  (I) 

NOQ  ==  ORDER  OF  GAUSSIAN  QUADRATURE  (I) 

PWT  ==  ARRAY  CONTAINING  PSI  AND  WEIGHTS  FOR  QUADRATURE  (I) 

PGF.NL  ==  VECTOR  OF  GLOBAL  EQUIVALENT  NODAL  LOADS  (0) 

PEENL  ==  EQUIVALENT  NODAL  LOADS  IN  ELEMENT  COORDINATES 
DN  ==  SHAPE  FUNCTION  MATRIX 
DL  ==  LENGTH  OF  ELEMENT 

P  ==  LOAD  VECTOR  FOR  INTERNAL  USE,  ELEMENT  COORDINATES 

CLEAR  PEENL  AND  ESTABLISH  P 

DO  10  1=1,6 

PEENL (1 5=0 . OD+OO 
10  CONTINUE 

P(1)=0. OD+OO 
P(2)=PR 
P(3)=0. OD+OO 
C 

0  PERFORM  GAUSSIAN  QUADRATURE  TO  GENERATE  EQUIVALENT  NODAL  LOADS 
C 

DO  100  1=1, NOQ 

CALL  NMATUl  ,R1 .  X2  ,R2  ,PWT(  1 , 1)  ,DN,DL) 

DO  90  J 1  =  1 ,6 
DO  90  J2=l,3 

PEENL (J 1 )=PEENL( J1)+(PWT(2 , I)*DN( J2, J1)*P( J2) ) 

90  CONTINUE 

100  CONTINUE 

DO  110  1=1,6 

PEENL  ( I )  =PEF.NL  ( I )  *DL/2 .  OD+OO 
110  CONTINUE 
C 

C  ROTATE  TO  GLOBAL  COORDINATES 
C 

CALL  ROTATE (X ' ,R1 , X2 , R2,6 , I , PEENL.2.PGENL) 

C 

C  FINISHED 
C 

RETURN 

END 

C************************************************************************** 

C  THIS  SUBROUTINE  INPUTS  GEOMETRY  AND  MATERIAL  PROPERTIES  AND  RETURNS  THE 
C  STRESS  RECOVERY  MATRIX.  THIS  MATRIX  WHEN  PREMULTIPLYING  GLOBAL 
C  DISPLACEMENTS  YIELDS  STRESSES  AT  THE  MIDPLANE  AND  +-  T/2  IN  ELEMENT 
C  COORDINATES  AT  THE  CENTER  OF  THE  ELEMENT. 

0 

SUBROUTINE  RECOVER  IX  1  ,K1  ,  X2 ,  R2  .  DNII,  E  ,T .  SRKOj 

DOUBLE  PRECISION  X  I  ,X2,K1  ,  R2  ,  DHU ,  E ,  T ,  SRI'.C  ( V  ,  t. )  ,G(5,5)  ,H(5,G)  , 

X  LI. ,  BMC  (7,5)  ,'I'INV.C  ;  .TEMP  (7 ,6) 


c 

C  EXPLANATION  OF  VARIABLES,  (I)  INPUT,  (0)  OUTPUT 
C  Xi,  Ri  ==  COORDINATES  OF  END  POINTS  (I) 

C  DNU  ==  POISSONS  RATIO  (I) 

C  F.  ==  MODULUS  (I) 

C  T  ==  ELEMENT  THICKNESS  (I) 

C  SREC  ==  STRESS  RECOVERY  MATRIX  (0) 

C  G  ==  MATERIAL  PROPERTY  MATRIX 

C  B  ==  B  MATRIX 

C  DL  ==  LENGTH  OF  ELEMENT 

C  BMC  ==  BENDING  MOMENT  CONVERSION  MATRIX 

C  TINV  ==  INVERSE  OF  THICKNESS 

C  Ci,  TEMPi  ==  DUMMY  CONSTANTS  ft  MATRICES 

C 

C  CLEAR  SREC,  TEMP  AND  BMC  MATRICES 
C 

DO  10  1=1,7 
DO  8  J=1 ,6 

TEMP (I , J) =0 .OD+OO 

8  CONTINUE 
DO  9  J=1 , 5 

BMC(I,J)=0. OD+OO 

9  CONTINUE 

10  CONTINUE 
C 

C  CREATE  CONSTANTS  FOR  BMC  CREATION 
C 

TINV=1 .OD+OO/T 
Cl=6 . OD+OO/ (T*T) 

C 

C  CREATE  BMC,  G,  AND  B  MATRICES 

C 

BMC( 1 , 1)=TINV 
BMC(1 ,3)=C1 
BMC(2,2)=TINV 
BMC(2 ,4)=C1 
BMC(3 , 1)=TINV 
BMC(4,2)=TINV 
BMC(5, 1)=TINV 
BMC(5,3)=-C1 
BMC(S,2)=TINV 
BMC(6,4)=-C1 
BMC(7,5)=1. OD+OO 

CALL  BMAT(X1,R1,X2,R2,0. OD+OO , B , DL) 

CALL  EMAT(DNU,E,T,G) 

C 

C  DO  BMC*E*B 
C 

DO  100  Jl  =  l  ,7 
DO  100  J2=l,5 
DO  100  J3=1,S 
DO  100  J4=l , 6 

TEMP(J1 , J4)=TEMP( J1 , J4)+BMC (J1,J2)*G(J2,J3)*B(J3,J4) 

100  CONTINUE 
C 

C  TRANSFORM  TO  ALLOW  OPERATION  OFF  GLOBAL  DISPLACEMENTS 
C 

CALL  ROTATE (X 1 , Rl , X2 ,R2 ,7 , 6 .TEMP , 3 .SREC) 

C 

C  FINISHED 
C 

RETURN 

END 

C** ************************************************************************ 

C  THIS  SUBROUTINE  INPUTS  GEOMETRY  AND  THE  GENERALIZED  COORDINATE  OF  INTEREST 
C  AND  RETURNS  THE  B  MATRIX  AND  ELEMENT  LENGTH. 

C 

SUBROUTINE  BMATCX 1 , Rl , X2 , R2 , DKSI ,B , DL) 

DOUBLE  PRECISION  XI  ,R1  ,X2,R2,DKSI,B(S,f>)  ,DL,RO ,PHI .CPR.SPR, R, 
ftDKM , DKP , DLINV 
C 

C  EXPLANATION  OF  VARIABLES,  (I)  INPUT,  (0)  OUTPUT 
C  X],  Ri  ==  COORDINATES  OF  END  POINTS  (I) 

C  DKSI  ==  KSI  TRANSFORMED  VARIABLE,  S  -->  KSI  C) 

C  B  ==  B  MATRIX,  DERIVATIVE  OF  SHAPE  FCTN ,  SO  SPEAK  (0) 

C  DL  ==  LENGTH  OF  ELEMENT  (0) 

C  RO  -  RADIUS  AT  CENTER  OF  ELEMENT 

C  PHI  =-  ANGLE  FROM  X,  AXIS  PER  COOK,  FIG  12.4-'.. 


C  CLEAR  D  MATRIX 
C 

DO  10  1  =  1  ,S 
DO  8  J=  1 ,6 

B (I , J) =0 . CD+OO 
8  CONTINUE 

10  CONTINUE 
C 

C  CALCULATE  RO,L,PHI 
C 

R0= (R2+R1) /2 . OD+OO 
PHI=DATAN2 (R2-R1 ,X2-X1) 

DL= (X2-X1) /DCOS(PHI) 

C 

C  CREATE  CONSTANTS  FOR  ELEMENT  CREATION 
C 

R=RO+(DKSI*DL*DSIN(PHI) / 2 . OD+OO) 

CPR=DCOS(PHI)/R 
SPR=DSIN(PHI)/R 
DKM=0 . 5D+00» ( 1 . OD+OO-DKSI) 

DKP=0 . 5D+00* ( 1 . OD+OO+DKSI) 

DLINV=1 . OD+OO/DL 
C 

C  NOW  CREATE  B  MATRIX 
C 

B(1 ,2)=-DLINV 
B ( 1 ,5)=DLINV 
B(2 ,2)=DKM*SPR 
B(2,5)=DKP*SPR 
B(2 , 1)=DKM»CPR 
B(2 ,4)=DKP*CPR 

B (2 , 3) =DL* ( 1 . OD+OO- (DKSI+DKSI) ) +CPR/8 . OD+OO 
B(2 , 6)=-B(2, 3) 

B(3,3)=DLINV 
B(3,6)=-DLINV 
B (4 ,  3)=-DKM*SPR 
B(4,6)=-DKP*SPR 
B(5, 1)=-DLINV 
B (5 ,4)=DLINV 

B (5 , 3) = (-DKSI/2 . OD+OO) -DK M 
B(5, 6)  =  (DKSI/2. OD+OO) -DKP 
C 

C  FINISHED 
C 

RETURN 

END 

C 

Cy*  *************************  ************************************************ 

C  THIS  SUBROUTINE  INPUTS  MATERIAL  PROPERTIES  AND  RETURNS  THE  MATERIAL 
C  MATRIX  FOR  USE  IN  CREATING  THE  STIFFNESS  MATRIX. 

C 

SUBROUTINE  EMAT(DNU,E,T,G) 

DOUBLE  PRECISION  DNU.E.T ,G (5 ,5) ,C1 
C 

C  EXPLANATION  OF  VARIABLES,  (I)  INPUT,  (0)  OUTPUT 
C  DNU  ==  POISSONS  RATIO  (I) 

C  E  ==  MODULUS  (I) 

C  T  ==  ELEMENT  THICKNESS  (I) 

C  G  ==  MATERIAL  PROPERTY  MATRIX  (0) 

r 

C  CLEAR  G  MATRIX 
C 

DO  10  1=1,5 
DO  8  J=1,S 

G(I, J)=0. OD+OO 
8  CONTINUE 

10  CONTINUE 

C 

C  CREATE  CONSTANTS  F^R  ELEMENT  CREATION 
C 

C 1=E»T/ ( 1 . OD+OO- (DNU+DNU) ) 

C 

C  NOW  CREATE  F.MAT  MATRIX 
C 

0(1  , 1)=C1 
0(1  , 2)  =  DNU+C 1 
0(2  ,  1  )=0( 1 ,2) 

0(2,2)=0] 

G(  '■  .3)=T«T*C1/12.(mm0ii 


oonooon  ooooo  non  ooo  onn^cc  oooonooooo  o  o  no  ooo 


G (3 ,4)=DNU*G(3 , 3) 

G(4,3)=G(3,4) 

G (4 ,4)=G(3 ,3) 

G(S,S)=5.0D+00*E*T/(12. 0D+00*( 1 . 0D+00+DNU)) 

FINISHED 

RETURN 

END 

************************************************************************** 

THIS  SUBROUTINE  INPUTS  GEOMETRY  AND  THE  GENERALIZED  COORDINATE  OF  INTEREST 
AND  RETURNS  THE  N  MATRIX  AND  ELEMENT  LENGTH 

SUBROUTINE  NMAT(X1 , R1 ,X2 ,R2 ,DKSI , DN , DL) 

DOUBIE  PRECISION  XI ,R1 ,X2 ,R2 ,DKSI ,DN (3 ,6) ,DL,DN1 ,DN2 ,DN3 

EXPLANATION  OF  VARIABLES,  (I)  INPUT,  (0)  OUTPUT 
Xi,  Ri  ==  COORDINATES  OF  END  POINTS  (I) 

DKSI  ==  KSI  TRANSFORMED  VARIABLE,  S  -->  KSI  (I) 

DN  ==  N,  THE  SHAPE  FUNCTION  MATRIX  (0) 

DL  ==  LENGTH  OF  ELEMENT  (0) 

DNi  ==  ELEMENTS  OF  SHAPE  FUNCTION  MATRIX 

CLEAR  DN  MATRIX 

DO  10  1=1,3 
DO  8  J=1 , 6 

DN(I , J)=0 .0D+00 
CONTINUE 
0  CONTINUE 

CALCULATE  DL 

DL=(X2-X1)/DC0S (DATAN2(R2-R1 ,X2-X1) ) 

DN1=0 . 5D+00* ( 1 . OD+OO-DKSI) 

DN2=0 . 5D+00* ( 1 . 0D+00+DKSI) 

DN3=DL* ( 1 . 0D+00-(DKSI*DKSI) ) /8 . 0D+00 

NOW  CREATE  DN  MATRIX 

DN( 1 ,2)=DN 1 
DN(1 ,S)=DN2 
DN(2 , 1)=DN1 
DN(2 ,4)=DN2 
DN(2,3)=DN3 
DN(2 ,6) =-DN3 
DN(3,3)=DN1 
DN(3 ,6)=DN2 

FINISHED 

RETURN 

END 

******************************************** ******************** ********** 

THIS  SUBROUTINE  INPUTS  GEOMETRY  AND  A  MATRIX  ENTITY,  CREATES  ROTATION 
MATRIX  USING  THE  FORM  u=TU  (CAPS  GLOBAL)  ,  AND  RETURNS  A  ROTATED  ENTITY  PER 
INSTRUCTION  CODE.  ENTITY  MUST  BE  COMPATIBLE  WITH  THE  OPERATION. 

SUBROUTINE  ROTATE (XI ,R1 ,X2 ,R2,NR0W,NC0L ,ENTIT, INSTR.ENOUT) 

DOUBLE  PRECISION  XI ,R1 , X2 ,R2  ,ENTIT(NROW,NCOL) ,T(6 ,6) ,PHI ,DC ,DS , 
ftENOUT(NROW , NCOL) 

EXPLANATION  OF  VARIABLES,  (I)  INPUT,  (0)  OUTPUT 
Xi,  Ri  ==  COORDINATES  OF  END  POINTS  (I) 

NROW  ==  NUMBER  OF  ROWS  IN  ENTIT  (I) 

NCOL  ==  NUMBER  OF  COLUMNS  IN  ENTIT  (I) 

ENTIT  ==  ENTITY  THAT  IS  BEING  TRANSFORMED  (I) 

INSTR  ==  INSTRUCTIONS  FOR  ROTATION  AS  BELOW  (I) 

C  1  — >  PREMULTIPLY  BY  T 

C  2  -->  PREMULTIPLY  BY  T  TRANSPOSE 

C  3  -->  POSTMULTIPLY  ENTIT  BY  T 

C  4  — >  POSTMULTIPLY  ENTIT  BY  T  TRANSPOSE 

C  S  -->  PRE  ft  POST  BY  T  TRANSPOSE  &  T,  RESPECTIVELY 

C  d  -->  PRE  a  POST  BY  T  a  T  TRANSPOSE,  REST ECTVIE1.Y 

C  F.NOUT  =*  TRANSFORMED  ENTITY  (0) 

C  T  ==  TRANSFORMATION  MATRIX  TO  GIVE  u=TU 

C  PHI  ==  ANGLE  FROM  AXIS  PER  COOK.  FIG  12.4-  3  » 

C 

C  CHECK  FOR  COMPATIBILITY 


;■  i 


<: 


I F(NCOL. NE.  6 . AND. NKOU . NE . 6)  CO  TO  510 
C 

C  CALCULATE  PHI,  COS(PHI),  SIN(PIII) 

C 

PHI=DATAN2(R2-R1,X2-X1) 

DC=DCOS(PHI) 

DS=DSIN (PHI)  < 

C  CLEAR  ROTATION  AND  ENTITY  OUTPUT  MATRICES 
C 

DO  10  1=1,6 
DO  6  .1=1,6 

T(I ,  J) =0 . OD+OO 
6  CONTINUE 

10  CONTINUE 

DO  20  1=1 ,NCOL 
DO  18  J=1 ,NROW 

F.N0UT(J,I)=0. OD+OO 
18  CONTINUE 

20  continue: 

C 

C  CREATE  ROTATION  MATRIX 
C 

T ( 1 , 1)=DC 
T(1,2)=-DS 
T(2, 1)=DS 
T (2 , 2)=DC 
T(4 ,4)=DC 
T(4,5)=-DS 
T(5 ,4) =DS 
T (5 ,5)=DC 
T(3 , 3)  =  1 . OD+OO 
T(6,6)=l .OD+OO 
C 

C  TRANSFORM  ENTITY 
C 

IF(INSTR.EQ.l)  THEN 
DO  100  J  1  =  1 , 6 
DO  100  J2=l ,6 

00  100  J3=l , NCOL 

ENOUTUl ,  J3)=EN0UT(I1 ,  J3)  +  (T(J1,J2)  *ENTIT (J2 ,  J3) ) 

100  CONTINUE 

ELSEIFCINSTR.EQ .2)  THEN 
DO  110  J 1= 1 ,6 
DO  110  J2=l,6 

DO  110  J3= 1 , NCOL 

EN0UT(J1.J3)=EN0UT(J1,J3)+(T(J2,J1)*ENTIT(J2,J3)) 

110  CONTINUE 

ELSEIFCINSTR.EQ .3)  THEN 
DO  120  J 1= 1 ,NROW 
DO  120  J2=! ,6 
DO  120  J3=l  ,6 

ENOUT( J1 , J3)=ENUUT( J1 , J3)+(ENTIT( J1 ,J3)*T(J2,J3)) 

120  CONTINUE 

ELSEIF(INSTR. EQ .4)  THEN 
DO  130  J 1 = 1 ,  NROW 
DO  130  J2= 1 ,6 
DO  130  J3=l ,6 

ENOUT  ( J 1 ,  J3)=EN0,n(  J 1 ,  J3)  +  (ENTIT(  J 1  ,J3)*T(J3,J2)) 

130  continue: 

ELSEIFCINSTR. EQ.S)  THEN 

IFCNCOL . NE . 6 . OR . NROW . NE . 6)  GO  TO  500 
DO  140  Jl  =  1  ,6 
DO  140  J2= 1 ,6 
DO  140  J3=1,G 
DO  140  .74  =  1  ,6 

ENOUT  (.5 1  ,  ,'4)=ENGUT  (.71 ,  ,74)+T(.72 ,  J1 7  »ENTIT(  J2,  J3)  *T(  J3  ,  J4) 

140  continue: 

EI.SE I  FC7 NSTR .  EQ  . 6 )  THEN 

I F (NCOL . NE . 6 . OR . NROW . NE . 6)  GO  TO  500 
DO  150  .71  =  1.6 
DO  160  .72=1  .6 
DO  150  ]3=1, 6 

DfJ  160  14=1,6  ,  ,  . 

ENOUT  7  7!  .  ! 17  ENOUT C  II  ,  14 )  +T(  II ,  72 ; -ENTITC  J2,  73 )  •  !  <  1 1  ,  4  '■) 

160  CONTINUE 


ELSE 

PKI 


! N  ROTATE,  incorrect  SIZING 


STOP 
ENDIF 
GO  TO  GOO 
C 

C  ERROR  STATEMENTS 
C 

500  PRINT*  , ’ERROR  IN  ROTATE,  F.NTIT  NOT  SQUARE’ 

STOP 

510  PRINT* , ’ERROR  IN  ROTATE,  ROWS  AND/OR  COLUMNS  NOT  COMPATIBLE’ 
STOP 


C 

C  FINISHED 
C 

600  RETURN 
END 

c*********$**************************************************************** 

C  THIS  SUBROUTINE  INFUTS  FILE  DATA  AND  WRITES  AN  ASTROS  COMPATIBLE  DIRECT 
C  MATRIX  INPUT  BULK  DATA  FILE  USING  LARGE  FIELDS. 

C 

SUBROUTINE  ASTOUT (MATRIX .PRECIS .FORM .NROWS, NCOLS , INFILE, NOFILE) 

DOUBLE  PRECISION  a(4) 

DIMENSION  INSTR(S)  ,  ICRI4) 

CHARACTER  MATRIX*8,  PRECIS*8,  F0RM*8,  INFILE*8 
OPEN (UNIT=8 ,  FILE=INFILE,  STATUS= ’OLD ’ ) 

REWIND(8) 

C 

C  EXPLANATION  OF  VARIABLES,  (I)  INPUT,  (0)  OUTPUT 
C  MATRIX  ==  CHARACTER  DENOTING  MATRIX  NAME  (I) 

C  PRECIS  ==  CHARACTER  DENOTING  PRECISION  (I) 

C  FORM  ==  CHARACTER  DENOTING  FORM  OF  MATRIX  (I) 

C  NROWS,  NCOLS  ==  NUMBER  OF  ROWS  AND  COLUMNS  IN  MATRIX  (I) 

C  INFILE  ==  INPUT  FILE  NAME  (I) 

C  NOFILE  ==  FORTRAN  OUTPUT  FILE  NUMBER  (I) 

C  A  ==  REAL  MATRIX  TERM 

C  INSTR  ==  INSTRUCTION  CODE  READ  FROM  INPUT  FILE  NUMBER  NIFILE 
C  1  — >  FOLLOWING  NUMBER  IS  REAL,  DOUBLE  PRECISION 

C  2  — >  FOLLOWING  NUMBER  IN  FILE  IS  INTEGER 

C  3  — >  THIS  IS  THE  END  OF  THE  FILE 

C  ICR  ==  INTEGER  COLUMN  OR  ROW  NUMBER  READ  FROM  INPUT  FILE  NIFILE 

C  IC2  ==  INSTRUCTION  CODE,  NUMBER  OF  VALUES  TO  OUTPUT,  1  — >  4 

C  IC3  ==  INSTRUCTION  CODE 

C  1  — >  THIS  IS  LAST  LINE  TO  OUTPUT 

C  2  -->  DON’T  STOP,  MORE  LINES  WILL  FOLLOW 

C 

C  ESTABLISH  POSSIBLE  FORMATS 
C 


1 

FORMAT! ’DMI 

’ , 3A8 ,218 ,24X , ’ABC’) 

2 

FORMAT! ’*BC 

’ , 4E16 . 8 , A3) 

3 

FORMAT !’*BC 

’ , I 16 , 3E16 . 8 , A3) 

4 

FORMAT! ’*BC 

’ , 21 16 ,2E16. 8 , A3) 

5 

FORMAT! ’*BC 

’ ,E16.8,2I16,E16.8,A3) 

6 

FORMAT !’*BC 

’ , 2E16 . 8 , 2116 , A3) 

7 

FORMAT !’*BC 

’ , 3E16 . 8 , 116  ,  A3) 

8 

FORMAT! ’*BC 

’ .4E16.8) 

9 

FORMAT! ’*BC 

’ .3E16.8) 

10 

FORMAT !’*BC 

’ .2E16.8) 

11 

FORMAT! ’*BC 

’ , E16 . 8) 

12 

FORMAT!  ’*BC 

’  ,  1 16 , 3E16 . 8) 

13 

FORMAT! ’*BC 

’  ,I16,2E16.8) 

14 

FORMAT! ’*BC 

’ , I 16 ,E16 . 8) 

15 

FORMAT!  ’*BC 

’ , 21 16 , 2E16 . 8) 

16 

FORMAT! ’*BC 

’ ,2116, E16. 8) 

17 

FORMAT! ’*BC 

’ ,E1G.8,2I16,E16.8) 

30 

FORMAT! 18) 

31 

F0RMAT!E16.8) 

C 

C  WRITE  FIRST  LINE 
C 

WRITE(NOFILF. ,  1)  MATRIX  .PRECIS  .FORM  .NROWS ,  NCOLS 
C 

C  STEP  THROUGH  NIFILE  AND  OUTPUT  PER  APPROPRIATE  FORMAT 


C 

100  READ (8 , 30)  INSTR( 1 ) 

IF(INSTR(1) .EQ. 3)  THEN 

PRINT* , ’FIRST  RECORD  READ  IN  RUB  ASTOUT  I?  END  11 
STOP 
ENDIF 

ISO  DO  200  r=l,4 

1 F ( INSTR ( I ) -Eg. 1)  THEN 


READ(fi ,  3  !  )  A  (  J  ) 

ELSEIF ( I NSTHC 1 )  .HQ. 2)  THEN 
RF.ADCo ,  30)  J  CR(l) 

ENDIF 

RFAD (8,30)  INSTRU+1) 

I F ( I NSTRCI  + 1 )  .F.Q.  3)  TUFN 
IC2=I 
IC3=1 
GO  TG  GOO 
KNDIF 

200  CONTINUE 
IC2=4 
IC3=2 

500  I F(IC2 . EQ . 1 )  THEN 

WRITE (NOFI !.F ,  11)  A(l) 

GO  TO  700 

EI.SEIF(IC2  .  Ft)  .2)  THEN 

IF(INSTRU)  .EQ.  1 .  AND.  INSTR(2)  .EQ.  1)  THEN 
WRITE (NOFILE, 10)  A(1),A(2) 

GO  TO  700 

ELSEIF(INSTRd)  .  EQ .  2  .  AND .  INSTR(2)  .  EQ  - 1 )  THEN 
WRITE(NOFILE.H)  ICR(1),A(2) 

GO  TO  700 
ENDIF 
GO  TO  600 

ELSEIFQC2.F.Q.3)  THEN 

IF(INSTRO)  .EQ.  1  .  AND .  INSTR(2)  . EQ .  1 .  AND.  INf.TR(3)  . EQ  .  1)  THEN 
WRITE (NO FILE ,9)  A ( 1 ) , A (2) , A (3) 

GO  TO  700 

ELSEIF(INSTR(1)  .  EQ  .  2  .  AND  .  INSTR(2)  .  EQ .  1 .  AND .  INSTR(3)  .EQ.  1)  THEN 
WRITE(NOFILE, 13)  ICR( 1 ) , A (2) , A (3) 

GO  TO  700 

ELSEIF(INSTR(1)  .F.Q  .  2 .  AND.  INSTR(2)  .  EQ .  2.  AND .  INSTRC3)  .EQ  .  1)  THEN 
WR;TE(N0FILE, 16)  ICR(l) ,ICR(2) ,A(3) 

GO  TO  700 
ENDIF 
GO  TO  601 

F.I.SEIF (IC2  .EQ  .4)  THEN 
IF(IC3 . EQ . 1 )  THEN 

IF(INSTR(1) . EQ . 1 . AND. INSTR(2) .EQ . 1 . AND. INSTR(3) .EQ. 1 .AND. 

&  INSTRd)  .  EQ  .  1 )  THEN 

WRITE(N0FILE,3)  A(l) ,A(2) ,A(3) ,A(4) 

CO  TO  700 

ELSEIF  (INS7R(1 )  .F.Q .  2 .  AND .  INSTRC2)  .EQ.  1 .  AND .  INSTR(3)  .EQ.  1  .AND. 
&  TNSTR(4) .Eq. 1)  THEN 

WRITE (NOFILE, 12)  ICR(1 ) , A (2) , A (3) , A (4 ) 

GO  TO  700 

EL  EIF(  INSTRd)  .  EQ .  2  .  AND .  INSTR(2)  .  EQ .  2  .  AND .  INSTR(3)  .EQ.  1 .  AND. 
&  INSTRd)  . EQ  .  1 )  THEN 

WRITE (NOFILE , 15)  ICR(l) ,ICR(2) ,A(3) ,A(4) 

GO  TO  700 

ELSEIF  (INSTRd  )  .  EQ .  1 .  AND  .  INSTRC2)  .  F,Q .  2 .  AND .  INSTR(3)  .  EQ .  2  .  AND  . 
&  INSTRd)  .F.Q.l)  THEN 

WRITE (NOFILE, 17)  A (1 ) , ICR(2) , ICR(3) .A (4) 

GO  TO  700 
ENDIF 
GO  TO  602 

KLSF.Ir  (IC3.E0.2)  THEN 

I  F(  INSTRd)  .EQ.  1  .  AND.  INSTR(2)  .EQ.  1 .  AND.  INSTRd)  .  F.Q  .  1  .AND. 
fc  INSTR(4) . Eq . 1 )  THEN 

WRITE  ( NL'F  I  i.E  ,2)  A(l)  ,A(2)  ,A(3)  ,A(4)  ,  ’AFC’ 

J*'STR(  1 )  =  I  NSTR(  5) 

GO  TO  150 

EI.SEI  F(  INSTRd)  .  EQ  .  2 .  AND  .  II  STR(2)  .F.Q.  1  .  AND.  i  NSTR(3)  .  EQ.  1  .  AND  . 
h  TNSTKd)  .EQ.  1)  THEN 

WRITE  I  NOE  I  I.E,  3)  iOR(l)  ,A(2)  ,A(3)  ,A(4 )  ,  'ARC  ’ 

INSTRCd-  INSTR(5) 

GO  TO  l‘d 

ELSE  I  F  ( I  KSTR(  I )  .F.Q.  2.  AND  I  NSTR(2)  .  KQ .  2 .  *•  ND .  I  I.'STR '  3 )  .F.Q.  !  .AN!'. 

n.  ;n::tr(4)  . eq.  :  >  then 

Wk f TK( NQF  ILK  ,4 )  !CK(  1 )  ,  I  GR(2)  ,  ACi)  ,  Ad  .  ,  ARC  ' 

J  N.iTR  I  1 )  -  I  ND'  i  [if  r.) 

GO  T  )  i 

ELSE. !  F  d  N::TRf  !  >  .  KQ  .  1  .  AND  .  1  HGT!:(2J  .  FQ  .  2  .  AN:.-.  J  NDTR '  3  )  .  FQ  .  2  AND 
k  I  NSTk'4  )  .  FQ  I)  TIIKN 

Wk  T!  d':l  ILL.1.)  A(l)  ,  I  Git  (2)  ,  ICI'.f  3)  ,A  ■  .  ,  Ai:0  • 
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&  1NSTIK4)  .EQ.2)  THEN 

WHITE (N0F1 I.E ,6)  A  C 1 ) ,  A  (2) ,  I CR(3) ,ICR(T) . ’ABC’ 

I NSTR( 1 ) =1 NSTR(S) 

GO  TO  ISO 

HI. SEIF (INSTRC  1 )  .  EQ .  1 .  AND .  INSTRC2) .  EQ .  1 . AND . I NSTR(3) . EQ . 1 . AND . 

&  INSTRU)  .EQ.2)  THEN 

WRITE ( NOE  1 1.E , 7 )  A(l)  , A (2)  ,A(3)  ,ICR(4)  ,  ’ARC" 

1NSTR( 1 )=1 NSTR (0) 

GO  TO  ICO 
ENUIF 
GO  TO  603 
ENDIF 
ENDIF 


ERROR  OUTPUT 
C 


600 

PRINT*, 

STOP 

’ERROR  WRITING 

601 

PRINT*, 

STOP 

’ERROR  WRITING 

602 

PRINT*, 

STOP 

’ERROR  WRITING 

603 

C 

PRINT*. 

STOP 

’ERROR  WRITING 

C  FINISHED 

C 

700 

RETURN 

END 

LAST  TWO  ELEMENT  LINE  IN  ASTOUT’ 

LAST  THREE  ELEMENT  LINE  IN  ASTOUT’ 
LAST  FOUR  ELEMENT  LINE  IN  ASTOUT’ 
INTERIM  FOUR  ELEMENT  LINE  IN  ASTOUT’ 
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